Nothing
#' @title Run a Genetic Algorithm to optimize a wind farm layout
#' @name genetic_algorithm
#' @description Run a Genetic Algorithm to optimize the layout of wind turbines
#' on a given area. The algorithm works with a fixed amount of turbines, a
#' fixed rotor radius and a mean wind speed value for every incoming wind
#' direction.
#'
#' @export
#'
#' @param Polygon1 The considered area as SpatialPolygon, SimpleFeature Polygon
#' or coordinates as matrix/data.frame
#' @param GridMethod Should the polygon be divided into rectangular or hexagonal
#' grid cells? The default is "Rectangular" grid. Hexagonal grids
#' are computed when assigning "h" or "hexagon" to this input variable.
#' @param Rotor The rotor radius in meter
#' @param n The amount of turbines
#' @param fcrR A numeric value used for grid spacing. Default is \code{5}
#' @param referenceHeight The height at which the incoming wind speeds were
#' measured. Default is \code{RotorHeight}
#' @param RotorHeight The height of the turbine hub
#' @param SurfaceRoughness A surface roughness length in meters.
#' With the terrain effect model, a surface roughness is calculated for every
#' grid cell using the elevation and land cover data. Default is \code{0.3}
#' @param sourceCCL The path to the Corine Land Cover raster (.tif). Only
#' required when the terrain effect model is activated.
#' @param sourceCCLRoughness The source to the adapted Corine Land Cover legend
#' as .csv file. Only required when terrain effect model is activated. As
#' default a .csv file within this package (\file{~/extdata}) is taken that
#' was already adapted manually.
#' @param Proportionality A numeric value used for the grid calculation, as it
#' determines the percentage a grid cell must overlay the area.
#' Default is \code{1}
#' @param iteration The number of iterations. Default is \code{20}
#' @param mutr A numeric mutation rate. Default is \code{0.008}
#' @param vdirspe A data.frame containing the wind speeds, directions and
#' probabilities. See \code{\link{windata_format}}.
#' @param topograp Boolean value, which indicates if the terrain effect model
#' should be enabled or not. Default is \code{FALSE}
#' @param elitism Boolean value, which indicates whether elitism should be
#' activated or not. Default is \code{TRUE}
#' @param nelit If \code{elitism} is TRUE, this input determines the amount
#' of individuals in the elite group. Default is 7
#' @param selstate Determines which selection method is used, "FIX" selects a
#' constant percentage and "VAR" selects a variable percentage, depending on
#' the development of the fitness values. Default is "FIX"
#' @param crossPart1 Determines which crossover method is used, "EQU" divides
#' the genetic code at equal intervals and "RAN" divides the genetic code at
#' random locations. Default is \code{"EQU"}
#' @param trimForce If \code{TRUE} the algorithm will use a probabilistic
#' approach to correct the windfarms to the desired amount of turbines.
#' If \code{FALSE} the adjustment will be random. Default is \code{FALSE}
#' @param Projection A spatial reference system. Depending on your PROJ-version,
#' it should either be a numeric `EPSG-code` or a `Proj4-string`.
#' Default is \code{EPSG:3035}
#' @param weibull A boolean value that specifies whether to take Weibull
#' parameters into account. If \code{TRUE}, the wind speed values
#' of \code{vdirspe} are ignored. The algorithm will calculate the mean
#' wind speed for every wind turbine according to the Weibull parameters.
#' Default is \code{FALSE}
#' @param weibullsrc A list of Weibull parameter rasters, where the first list
#' item must be the shape parameter raster `k` and the second item must be the
#' scale parameter raster `a` of the Weibull distribution. If no list is given,
#' then rasters included in the package are used instead, which currently only
#' cover Austria. This variable is only used if \code{weibull = TRUE}.
#' @param Parallel A boolean value, indicating whether parallel processing should
#' be used. The *parallel* and *doParallel* packages are used for parallel
#' processing. Default is \code{FALSE}
#' @param numCluster If \code{Parallel} is TRUE, this variable defines the
#' number of clusters to be used. Default is \code{2}
#' @param verbose If TRUE it will print information for every generation.
#' Default is \code{FALSE}
#' @param plotit If TRUE it will plot the best windfarm of every generation.
#' Default is \code{FALSE}
#'
#' @family Genetic Algorithm Functions
#' @return The result is a matrix with aggregated values per generation; the
#' best individual regarding energy and efficiency per generation, some fuzzy
#' control variables per generation, a list of all fitness values per
#' generation, the amount of individuals after each process, a matrix of all
#' energy, efficiency and fitness values per generation, the selection and
#' crossover parameters, a matrix with the generational difference in maximum
#' and mean energy output, a matrix with the given inputs, a dataframe with
#' the wind information, the mutation rate per generation and a matrix with all
#' tested wind farm layouts.
#'
#' @details A terrain effect model can be included in the optimization process.
#' Therefore, a digital elevation model will be downloaded automatically via the
#' \code{elevatr::get_elev_raster} function. A land cover raster can also be
#' downloaded automatically from the EEA-website, or the path to a raster file
#' can be passed to \code{sourceCCL}. The algorithm uses an adapted version of the
#' Raster legend ("clc_legend.csv"), which is stored in the package directory
#' \file{~/inst/extdata}. To use other values for the land cover roughness
#' lengths, insert a column named \strong{"Rauhigkeit_z"} to the .csv file,
#' assign a surface roughness length to all land cover types. Be sure that all
#' rows are filled with numeric values and save the file with \strong{";"}
#' separation. Assign the path of the file to the input variable
#' \code{sourceCCLRoughness} of this function.
#'
#' @examples \donttest{
#' ## Create a random rectangular shapefile
#' library(sf)
#'
#' Polygon1 <- sf::st_as_sf(sf::st_sfc(
#' sf::st_polygon(list(cbind(
#' c(4498482, 4498482, 4499991, 4499991, 4498482),
#' c(2668272, 2669343, 2669343, 2668272, 2668272)))),
#' crs = 3035
#' ))
#'
#' ## Create a uniform and unidirectional wind data.frame and plot the
#' ## resulting wind rose
#' data.in <- data.frame(ws = 12, wd = 0)
#' windrosePlot <- plot_windrose(data = data.in, spd = data.in$ws,
#' dir = data.in$wd, dirres=10, spdmax=20)
#'
#' ## Runs an optimization run for 20 iterations with the
#' ## given shapefile (Polygon1), the wind data.frame (data.in),
#' ## 12 turbines (n) with rotor radii of 30m (Rotor) and rotor height of 100m.
#' result <- genetic_algorithm(Polygon1 = Polygon1,
#' n = 12,
#' vdirspe = data.in,
#' Rotor = 30,
#' RotorHeight = 100)
#' plot_windfarmGA(result = result, Polygon1 = Polygon1)
#'}
genetic_algorithm <- function(Polygon1, GridMethod, Rotor, n, fcrR, referenceHeight,
RotorHeight, SurfaceRoughness, Proportionality,
iteration, mutr, vdirspe, topograp, elitism, nelit,
selstate, crossPart1, trimForce, Projection,
sourceCCL, sourceCCLRoughness, weibull, weibullsrc,
Parallel, numCluster, verbose = FALSE, plotit = FALSE){
## set Graphic Params ###############
if (plotit) {
oldpar <- graphics::par(no.readonly = TRUE)
on.exit(par(oldpar))
plot.new()
graphics::par(ask = FALSE)
}
## MISSING ARGUMENTS ###############
if (missing(fcrR)) {
fcrR <- 5
}
if (missing(topograp)) {
topograp <- FALSE
}
if (missing(GridMethod)) {
GridMethod <- "Rectangular"
}
if (missing(Parallel)) {
Parallel <- FALSE
}
if (missing(numCluster)) {
numCluster <- 2
}
if (missing(weibull)) {
weibull <- FALSE
}
if (missing(selstate)) {
selstate <- "FIX"
}
if (missing(crossPart1)) {
crossPart1 <- "EQU"
}
if (missing(SurfaceRoughness)) {
SurfaceRoughness <- 0.3
}
if (missing(Proportionality)) {
Proportionality <- 1
}
if (missing(mutr)) {
mutr <- 0.008
}
if (missing(elitism)) {
elitism <- TRUE
if (missing(nelit)) {
nelit <- 7
}
}
if (missing(trimForce)) {
trimForce <- FALSE
}
if (missing(RotorHeight)) {
stop("The variable 'RotorHeight' is not defined. Assign the turbine heights to 'RotorHeight'.")
}
if (missing(referenceHeight)) {
referenceHeight <- RotorHeight
}
if (missing(iteration)) {
iteration <- 20
}
if (missing(Projection)) {
if (utils::compareVersion(sf::sf_extSoftVersion()[[3]], "6") > 0) {
ProjLAEA <- 3035
} else {
ProjLAEA <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
}
} else {
ProjLAEA <- Projection
}
if (missing(vdirspe)) {
stop("No Winddata is given.")
}
if (missing(n)) {
stop("The variable 'n' is not defined. Assign the number of turbines to 'n'.")
}
if (missing(Rotor)) {
stop("The variable 'Rotor' is not defined. Assign the rotor radius to 'Rotor'.")
}
## INIT VARIABLES 1 #################
selstate <- toupper(selstate)
crossPart1 <- toupper(crossPart1)
## Is the Polygon Spatial / SF / coordinates - It will transform to sf-Polygons
Polygon1 <- isSpatial(Polygon1, ProjLAEA)
if (is.na(st_crs(Polygon1))) {
stop("The input area is not projected.")
}
## Grid size calculation
resol2 <- fcrR * Rotor
## Max Amount of individuals in the Crossover-Method
# CrossUpLimit <- 300
CrossUpLimit <- getOption("windfarmGA.max_population")
## Start Parallel Cluster ###############
## Is Parallel processing activated? Check the max number of cores and set to max-1 if value exceeds.
if (Parallel) {
max_cores <- parallel::detectCores()
if (numCluster > max_cores) {
warning("Maximum number of cores is: ", max_cores, "\n'numCluster' will be set to: ", max_cores-1)
numCluster <- max_cores - 1
}
type_cluster <- "PSOCK" ## TODO - should this be available as option too?
cl <- parallel::makeCluster(numCluster, type = type_cluster)
doParallel::registerDoParallel(cl)
}
## WEIBULL ###############
## Is Weibull activated? If no source is given, take values from package
if (weibull) {
if (verbose) {
cat("\nWeibull Distribution is used.")
}
if (missing(weibullsrc)) {
if (verbose) cat("\nWeibull Information from package is used.\n")
if (!file.exists("k120_100m_Lambert.tif") && !file.exists("a120_100m_Lambert.tif")) {
download.file("http://github.com/YsoSirius/windfarm_data/raw/master/weibulldata.zip",
destfile = "weibulldata.zip",
method = "auto")
unzip("weibulldata.zip")
unlink("weibulldata.zip")
}
k_weibull <- raster("k120_100m_Lambert.tif")
a_weibull <- raster("a120_100m_Lambert.tif")
} else {
if (verbose) {
cat("\nWeibull Information from input is used.\n")
}
## Project Shapefile to raster, Crop/Mask and project raster back
k_weibull <- weibullsrc[[1]]
a_weibull <- weibullsrc[[2]]
}
## Project Shapefile to raster proj, Crop/Mask and project raster back
shape_project <- st_transform(Polygon1, crs = st_crs(a_weibull))
k_par_crop <- raster::crop(x = k_weibull, y = raster::extent(shape_project))
a_par_crop <- raster::crop(x = a_weibull, y = raster::extent(shape_project))
weibl_k <- raster::mask(x = k_par_crop, mask = shape_project)
weibl_a <- raster::mask(x = a_par_crop, mask = shape_project)
estim_speed_raster <- weibl_a * (gamma(1 + (1 / weibl_k)))
estim_speed_raster <- raster::projectRaster(estim_speed_raster,
crs = raster::crs(Polygon1))
} else {
estim_speed_raster <- FALSE
}
## CHECK INPUTS ###############
## Check if Input Data is correct and prints it out.
if (crossPart1 != "EQU" & crossPart1 != "RAN") {
crossPart1 <- readinteger()
}
if (selstate != "FIX" & selstate != "VAR") {
selstate <- readintegerSel()
}
inputData <- list(
Input_Data = rbind("Rotorradius" = Rotor,
"Number of turbines" = n,
"Grid Shape Factor" = fcrR,
"Iterations" = iteration,
"Mutation Rate" = mutr,
"Percentage of Polygon" = Proportionality,
"Topographie" = topograp,
"Elitarism" = elitism,
"Selection Method" = selstate,
"Trim Force Method Used" = trimForce,
"Crossover Method Used" = crossPart1,
"Reference Height" = referenceHeight,
"Rotor Height" = RotorHeight,
"Resolution" = resol2,
"Parallel Processing" = Parallel,
"Number Clusters" = numCluster,
"Active Weibull" = weibull,
"Grid Method" = GridMethod,
"Projection" = ProjLAEA))
inputWind <- list(Windspeed_Data = vdirspe)
if (verbose) {
print(inputData)
print(inputWind)
}
## Winddata Formatting #######################
winddata <- windata_format(vdirspe)
#######################
## Project Polygon ###############
if (utils::compareVersion(sf::sf_extSoftVersion()[[3]], "6") > 0) {
if (suppressWarnings(!isTRUE( all.equal(st_crs(Polygon1), st_crs(ProjLAEA)) ))) {
Polygon1 <- sf::st_transform(Polygon1, ProjLAEA)
}
} else {
if (as.character(raster::crs(Polygon1)) != ProjLAEA) {
Polygon1 <- sf::st_transform(Polygon1, ProjLAEA)
}
}
## Make GRID ###############
## Calculate a Grid and an indexed data.frame with coordinates and grid cell Ids.
GridMethod <- toupper(GridMethod)
## Decide if the space division should be rectangular or in hexagons.
if (GridMethod != "HEXAGON" & GridMethod != "H") {
# Calculate a Grid and an indexed data.frame with coordinates and grid cell Ids.
Grid1 <- grid_area(Polygon1, resol2, Proportionality)
Grid <- Grid1[[1]]
grid_filtered <- Grid1[[2]]
} else {
# Calculate a Grid with hexagonal grid cells
Grid1 <- hexa_area(Polygon1, resol2)
Grid <- Grid1[[1]]
grid_filtered <- Grid1[[2]]
}
n_gridcells <- nrow(Grid)
## INIT VARIABLES 2 ###############
## Determine the amount of initial individuals and create initial population.
nStart <- (n_gridcells * n) / iteration
if (nStart < 100) {
nStart <- 100
}
if (nStart > 300) {
nStart <- 300
}
nStart <- ceiling(nStart)
startsel <- init_population(Grid, n, nStart)
## Initialize all needed variables as list.
maxParkwirkungsg <- 0
allparkcoeff <- vector("list", iteration)
bestPaEn <- vector("list", iteration)
bestPaEf <- vector("list", iteration)
fuzzycontr <- vector("list", iteration)
fitnessValues <- vector("list", iteration)
nindiv <- vector("list", iteration)
clouddata <- vector("list", iteration)
selcross <- vector("list", iteration)
beorwor <- vector("list", iteration)
mut_rate <- vector("list", iteration)
allCoords <- vector("list", iteration)
## TERRAIN EFFECT MODEL ###############
## Checks if terrain effect model is activated, and makes necessary caluclations.
if (!topograp) {
if (verbose) {
cat("Topography and orography are not taken into account.\n")
}
srtm_crop <- ""
cclRaster <- ""
} else {
if (verbose) {
cat("Topography and orography are taken into account.\n")
}
if (plotit) {
par(mfrow = c(3, 1))
}
if (missing(sourceCCL)) {
message("\nNo land cover raster ('sourceCCL') was given. It will be downloaded from ",
"the EEA-website.\n")
if (!file.exists("g100_06.tif")) {
# "https://www.eea.europa.eu/data-and-maps/data/clc-2006-raster-3/clc-2006-100m/g100_06.zip/at_download/file"
download.file("http://github.com/YsoSirius/windfarm_data/raw/master/clc.zip",
destfile = "clc.zip",
method = "auto")
unzip("clc.zip")
unlink("clc.zip")
}
ccl <- raster::raster("g100_06.tif")
} else {
ccl <- raster::raster(sourceCCL)
}
cclPoly <- raster::crop(ccl, Polygon1)
cclPoly1 <- raster::mask(cclPoly, Polygon1)
## SRTM Daten
Polygon_wgs84 <- sf::st_transform(Polygon1, st_crs(4326))
# extpol <- round(apply(matrix(st_bbox(Polygon1), ncol = 2), 1, mean))
# srtm <- tryCatch(raster::getData("SRTM", lon = extpol[1], lat = extpol[2]),
# error = function(e) {
# stop("\nCould not download SRTM for the given Polygon.",
# "Check the Projection of the Polygon.\n", call. = FALSE)
# })
srtm <- tryCatch(elevatr::get_elev_raster(verbose = verbose,
locations = as(Polygon_wgs84, "Spatial"), z = 11),
error = function(e) {
stop("\nDownloading Elevation data failed for the given Polygon.\n",
e, call. = FALSE)
})
srtm_crop <- raster::crop(srtm, Polygon_wgs84)
srtm_crop <- raster::mask(srtm_crop, Polygon_wgs84)
srtm_crop <- raster::projectRaster(srtm_crop, crs = raster::crs(Polygon1))
if (plotit) {
raster::plot(srtm_crop, main = "Elevation Data")
plot(Polygon1, add = TRUE, color = "transparent")
plot(grid_filtered, add = TRUE)
}
roughrast <- raster::terrain(srtm_crop, "roughness")
if (all(is.na(raster::values(roughrast)))) {
warning("Cannot calculate a surface roughness. \nMaybe the resolution or ",
"the area is too small. Roughness values are set to 1.\n")
raster::values(roughrast) <- 1
}
srtm_crop <- list(
strm_crop = srtm_crop,
orogr1 = raster::calc(srtm_crop, function(x) {
x / (raster::cellStats(srtm_crop, mean, na.rm = TRUE))
}),
roughness = roughrast
)
# Include Corine Land Cover Raster to get an estimation of Surface Roughness
if (missing(sourceCCLRoughness)) {
path <- paste0(system.file(package = "windfarmGA"), "/extdata/")
sourceCCLRoughness <- paste0(path, "clc_legend.csv")
} else {
if (verbose) {
message("You are using your own Corine Land Cover legend.\n")
}
}
rauhigkeitz <- utils::read.csv(sourceCCLRoughness,
header = TRUE, sep = ";")
cclRaster <- raster::reclassify(cclPoly1,
matrix(c(rauhigkeitz$GRID_CODE,
rauhigkeitz$Rauhigkeit_z),
ncol = 2))
if (plotit) {
raster::plot(cclRaster, main = "Surface Roughness from Corine Land Cover")
}
}
## GENETIC ALGORITHM #################
if (verbose) {cat("\nStart Genetic Algorithm ...\n")}
rbPal <- grDevices::colorRampPalette(c("red", "green"))
i <- 1
while (i <= iteration) {
if (!verbose) {
cat(".")
}
## FITNESS (+getRectV) ###############
if (i == 1) {
fit <- fitness(selection = startsel, referenceHeight = referenceHeight,
RotorHeight = RotorHeight,
SurfaceRoughness = SurfaceRoughness,
Polygon = Polygon1, resol1 = resol2, rot = Rotor,
dirspeed = winddata, srtm_crop = srtm_crop,
topograp = topograp, cclRaster = cclRaster,
weibull = estim_speed_raster,
Parallel = Parallel, numCluster = numCluster)
} else {
getRectV <- get_grids(mut1, Grid)
fit <- fitness(selection = getRectV, referenceHeight = referenceHeight,
RotorHeight = RotorHeight,
SurfaceRoughness = SurfaceRoughness,
Polygon = Polygon1, resol1 = resol2, rot = Rotor,
dirspeed = winddata, srtm_crop = srtm_crop,
topograp = topograp, cclRaster = cclRaster,
weibull = estim_speed_raster,
Parallel = Parallel, numCluster = numCluster)
}
## Fitness Result Processing ###############
allparks <- do.call("rbind", fit)
allparksUni <- subset.matrix(allparks,
subset = !duplicated(allparks[, "Run"]))
allCoords[[i]] <- allparks
maxparkfitness <- round(max(allparksUni[, "Parkfitness"]), 4)
meanparkfitness <- round(mean(allparksUni[, "Parkfitness"]), 3)
minparkfitness <- round(min(allparksUni[, "Parkfitness"]), 3)
MaxEnergyRedu <- round(max(allparksUni[, "EnergyOverall"]), 2)
MeanEnergyRedu <- round(mean(allparksUni[, "EnergyOverall"]), 2)
MinEnergyRedu <- round(min(allparksUni[, "EnergyOverall"]), 2)
maxParkwirkungsg <- round(max(allparksUni[, "EfficAllDir"]), 2)
meanParkwirkungsg <- round(mean(allparksUni[, "EfficAllDir"]), 2)
minParkwirkungsg <- round(min(allparksUni[, "EfficAllDir"]), 2)
allparkcoeff[[i]] <- cbind(
maxparkfitness, meanparkfitness, minparkfitness,
MaxEnergyRedu, MeanEnergyRedu, MinEnergyRedu,
maxParkwirkungsg, meanParkwirkungsg, minParkwirkungsg)
clouddata[[i]] <- subset.matrix(allparksUni,
select = c("EfficAllDir",
"EnergyOverall",
"Parkfitness"))
if (verbose) {
cat(c("\n\n", i, ": Round with coefficients ", allparkcoeff[[i]], "\n"))
}
## Highest Energy Output
xd <- max(allparks[, "EnergyOverall"])
ind <- allparks[, "EnergyOverall"] == xd
bestPaEn[[i]] <- allparks[ind, ][1:n, ]
## Highest Efficiency
xd1 <- max(allparks[, "EfficAllDir"])
ind1 <- allparks[, "EfficAllDir"] == xd1
bestPaEf[[i]] <- allparks[ind1, ][1:n, ]
# Print out most relevant information on Generation i
afvs <- allparks[allparks[, "EnergyOverall"] == max(
allparks[, "EnergyOverall"]), ]
if (verbose) {
cat(paste("How many individuals exist: ", length(fit) ), "\n")
cat(paste("How many parks are in local Optimum: ",
(length(afvs[, 1]) / n) ), "\n")
}
nindivfit <- length(fit)
if (plotit) {
lebre <- length(unique(bestPaEn[[i]][, "AbschGesamt"]))
if (lebre < 2) {
Col <- "green"
} else {
Col <- rbPal(lebre)[as.numeric(cut(-bestPaEn[[i]][, "AbschGesamt"],
breaks = lebre))]
}
lebre2 <- length(unique(bestPaEf[[i]][, "AbschGesamt"]))
if (lebre2 < 2) {
Col1 <- "green"
} else {
Col1 <- rbPal(lebre2)[as.numeric(cut(-bestPaEf[[i]][, "AbschGesamt"],
breaks = lebre2))]
}
}
x <- round(bestPaEn[[i]][, "EnergyOverall"][[1]], 2)
y <- round(bestPaEn[[i]][, "EfficAllDir"][[1]], 2)
e <- bestPaEn[[i]][, "EfficAllDir"]
x1 <- round(bestPaEf[[i]][, "EnergyOverall"][[1]], 2)
y1 <- round(bestPaEf[[i]][, "EfficAllDir"][[1]], 2)
e1 <- bestPaEf[[i]][, "EfficAllDir"]
##################
if (plotit) {
graphics::par(mfrow = c(1, 2))
plot(st_geometry(Polygon1), col="lightblue",
main = paste(i, "Round \n Best Energy Output: ", x,
"W/h \n Efficiency: ", y, "%"),
sub = paste("\n Number of turbines: ", length(e)))
plot(grid_filtered, add=TRUE)
graphics::points(bestPaEn[[i]][, "X"], bestPaEn[[i]][, "Y"],
col = Col, pch = 20, cex = 1.5)
plot(st_geometry(Polygon1), col="lightblue",
main = paste(i, "Round \n Best Efficiency Output: ",
x1, "W/h \n Efficiency: ", y1, "%"),
sub = paste("\n Number of turbines: ", length(e1)))
plot(grid_filtered, add = TRUE)
graphics::points(bestPaEf[[i]][, "X"], bestPaEf[[i]][, "Y"],
col = Col1, pch = 20, cex = 1.5)
}
## Fuzzy Control ###############
if (i > 20) {
besPE <- do.call("rbind", lapply(bestPaEn[1:i], function(x) {
max(x[, "EnergyOverall"])
}))
maxBisher <- max(besPE)
WhichMaxBs <- which(besPE == max(besPE))
if (length(WhichMaxBs) >= 2) {
BestForNo <- bestPaEn[sample(WhichMaxBs, 2)]
BestForNo[[1]][, "Run"] <- length(fit) + 1
BestForNo[[2]][, "Run"] <- length(fit) + 2
} else {
BestForNo <- bestPaEn[WhichMaxBs]
BestForNo <- append(BestForNo, BestForNo)
BestForNo[[1]][, "Run"] <- length(fit) + 1
BestForNo[[2]][, "Run"] <- length(fit) + 2
}
last7 <- besPE[i:(i - 5)]
if (!any(last7 == maxBisher)) {
if (verbose) {
cat(paste("Park with highest Fitness level to date ",
"is replaced in the list.", "\n\n"))
}
fit <- append(fit, BestForNo)
}
}
if (i == 1) {
## TODO I do have such a matrix already with that info or??
t0 <- subset.matrix(allparks, !duplicated(allparks[, "Run"]))
t0 <- t0[, "Parkfitness"]
fitnessValues[[i]] <- t0
rangeFitnessVt0 <- range(t0)
maxt0 <- max(t0)
meant0 <- mean(t0)
allcoef0 <- c(rangeFitnessVt0, meant0)
fuzzycontr[[i]] <- rbind(allcoef0)
colnames(fuzzycontr[[i]]) <- c("Min", "Max", "Mean")
teil <- 2
if (selstate == "VAR") {
teil <- 1.35
}
u <- 1.1
beorwor[[i]] <- cbind(0, 0)
}
## TODO - better fuzzy model and expose weights and values by options?
if (i >= 2 && i <= iteration) {
t0 <- subset.matrix(allparks, !duplicated(allparks[, "Run"]))
t0 <- t0[, "Parkfitness"]
fitnessValues[[i]] <- t0
rangeFitnessVt0 <- range(t0)
maxt0 <- max(t0)
meant0 <- mean(t0)
mint0 <- min(t0)
t1 <- fitnessValues[[i - 1]]
rangeFitnessVt1 <- range(t1)
maxt1 <- max(t1)
meant1 <- mean(t1)
mint1 <- min(t1)
maxDif <- maxt0 - maxt1
meanDif <- meant0 - meant1
minDif <- mint0 - mint1
WeightDif <- c(0.80, 0.2, 0.0)
maxunt <- (maxDif * WeightDif[1]) +
(meanDif * WeightDif[2]) + (minDif * WeightDif[3])
allcoef1 <- c(rangeFitnessVt0, meant0)
allcoef2 <- c(rangeFitnessVt1, meant1)
fuzzycontr[[i]] <- rbind(allcoef1, allcoef2)
colnames(fuzzycontr[[i]]) <- c("Min", "Max", "Mean")
if (maxunt < 0) {
pri <- "deteriorated"
teil <- teil - 0.02
u <- u - 0.06
} else if (maxunt == 0) {
pri <- "not changed"
teil <- teil; u <- u
} else {
pri <- "improved"
teil <- teil + 0.017
u <- u + 0.03
}
if (teil > 5) {
teil <- 5; u <- u + 0.09
if (verbose) cat("Min 20% Selected"); cat(paste("CPR is increased! CPR:", u, "SP: ", teil, "\n"))
}
if (trunc(u) < 0) {
u <- 0.5; teil <- teil - 0.4
if (verbose) cat(paste("Min 1 CrossPoints. Selection decreased. CPR:",u, "SP: ", teil, "\n"))
}
if (u >= 4) {
u <- 4; teil <- 4
if (verbose) cat(paste("Max 5 CrossPoints. Select fittest 25%. SP: ", teil, "\n"))
}
if (teil <= 4 / 3) {
teil <- 4 / 3
if (verbose) cat(paste("Max 75% selected. SP: ", teil, "\n"))
}
if (length(fit) <= 20) {
teil <- 1; u <- u + 0.1
if (verbose) cat(paste("Less than 20 individuals. Select all and increase ",
"Crossover-point rate. CPR: ", u, "SP: ", teil, "\n"))
}
# if (length(fit) <= 10) {
# teil <- 1; u <- u + 0.4
# if (verbose) cat(paste("Less than 10 individuals. Select all and increase ",
# "Crossover-point rate. CPR: ", u, "SP: ", teil, "\n"))
# }
# if (teil > 5) {
# teil <- 5
# if (verbose) cat(paste("Teil is bigger than 5. Set to max 5. SP:", teil, "\n"))
# }
u <- round(u, 2)
teil <- round(teil, 3)
if (verbose) {
cat(paste("Fitness of this population (", i,
"), compared to the prior,", pri,
"by", round(maxunt, 2), "W/h \n"))
}
meanunt <- meant0 - meant1
beorwor[[i]] <- cbind(maxunt, meanunt)
}
## SELECTION #################
if (selstate == "FIX") {
if (teil == 1) {
teil <- 1
} else {
teil <- 2
}
}
if (crossPart1 == "EQU") {
u <- round(u, 2)
}
## How many are selected and how much crossover points are used?
selcross[[i]] <- cbind(cross = trunc(u + 1), teil)
selec6best <- selection(fit = fit, Grid = Grid, teil = teil,
elitism = elitism, nelit = nelit,
selstate = selstate, verbose = verbose)
selec6best_bin <- selec6best[[1]]
if (verbose) {
cat(paste("Selection - Amount of Individuals: ",
length(selec6best_bin[1, -1]), "\n"))
}
nindivsel <- length(selec6best_bin[1, -1])
## CROSSOVER #################
## u determines the amount of crossover points,
## crossPart determines the method used (Equal/Random),
## uplimit is the maximum allowed permutations
crossOut <- crossover(se6 = selec6best, u = u, uplimit = CrossUpLimit,
crossPart = crossPart1,
verbose = verbose, seed = NULL)
if (verbose) {
cat(paste("Crossover - Amount of Individuals: ",
length(crossOut[1, ])))
}
nindivcros <- length(crossOut[1, ])
## MUTATION #################
## Variable Mutation Rate is activated if more than 2 individuals
## represent the current best solution.
loOp <- (length(afvs[, 1]) / n)
if (loOp > 2) {
mutrn <- round(runif(1, 0.03, 0.1), 2)
t1 <- (loOp * 1.25) / 42
mutrn <- mutrn * (1 + t1)
mutrn <- round(mutrn + ((i) / (20 * iteration)), 5)
mut <- mutation(a = crossOut, p = mutrn, seed = NULL)
mut_rat <- mutrn
if (verbose) {
cat(paste("\nVariable Mutation Rate is", mutrn, "\n"))
}
} else {
mut <- mutation(a = crossOut, p = mutr, seed = NULL)
mut_rat <- mutr
}
mut_rate[[i]] <- mut_rat
if (verbose) {
cat(paste("\nMutation - Amount of Individuals: ", length(mut[1, ])))
}
nindivmut <- length(mut[1, ])
## TRIMTON #################
## After Crossover and Mutation, the amount of turbines in a windpark
## change and have to be corrected to the required amount of turbines.
mut1 <- trimton(mut = mut, nturb = n, allparks = allparks,
nGrids = n_gridcells, trimForce = trimForce,
seed = NULL)
if (verbose) {
cat(paste("\nTrimToN - Amount of Individuals: ",
length(mut1[1, ])))
}
nindiv[[i]] <- cbind(nindivfit, nindivsel, nindivcros, nindivmut)
if (maxParkwirkungsg == 100) {
i <- iteration + 1
} else {
i <- i + 1
}
}
## Remove Parallel Cluster ###############
if (Parallel) {
parallel::stopCluster(cl)
# try(rm(cl), silent = TRUE)
}
## Reduce list, if algorithm didnt run all iterations. (Found Optimum) #################
mut_rate <- mut_rate[lapply(mut_rate, length) != 0]
beorwor <- beorwor[lapply(beorwor, length) != 0]
selcross <- selcross[lapply(selcross, length) != 0]
clouddata <- clouddata[lapply(clouddata, length) != 0]
allparkcoeff <- allparkcoeff[lapply(allparkcoeff, length) != 0]
bestPaEn <- bestPaEn[lapply(bestPaEn, length) != 0]
bestPaEf <- bestPaEf[lapply(bestPaEf, length) != 0]
fuzzycontr <- fuzzycontr[lapply(fuzzycontr, length) != 0]
fitnessValues <- fitnessValues[lapply(fitnessValues, length) != 0]
nindiv <- nindiv[lapply(nindiv, length) != 0]
allCoords <- allCoords[lapply(allCoords, length) != 0]
## Bind the results together and Output them. #################
alldata <- cbind(allparkcoeff, bestPaEn, bestPaEf,
fuzzycontr, fitnessValues, nindiv,
clouddata, selcross, beorwor,
inputData, inputWind, mut_rate, allCoords)
return(alldata)
}
#' @title Transform to Simple Feature Polygons
#' @name isSpatial
#' @description Helper Function, which transforms SpatialPolygons or coordinates
#' in matrix/data.frame - form to a Simple Feature Polygon
#'
#' @export
#'
#' @param shape An area as SpatialPolygon, SimpleFeature Polygon or coordinates
#' as matrix/data.frame
#' @param proj Which Projection should be assigned to matrix / data.frame
#' coordinates
#'
#' @family Helper Functions
#' @return A Simple Feature Polygon
#'
#' @details If the columns are named, it will look for common abbreviation to
#' match x/y or long/lat columns. If the columns are not named, the first 2
#' numeric columns are taken.
#'
#' @examples \donttest{
#' library(sf)
#' df <- rbind(c(4498482, 2668272), c(4498482, 2669343),
#' c(4499991, 2669343), c(4499991, 2668272))
#' isSpatial(df)
#'
#' Polygon1 <- sf::st_as_sf(sf::st_sfc(
#' sf::st_polygon(list(cbind(
#' c(4498482, 4498482, 4499991, 4499991, 4498482),
#' c(2668272, 2669343, 2669343, 2668272, 2668272)))),
#' crs = 3035
#' ))
#' isSpatial(st_coordinates(Polygon1), 3035)
#'}
isSpatial <- function(shape, proj) {
if (inherits(shape, "Spatial")) {
# shape <- as(shape, "Spatial")
shape <- st_as_sf(shape)
## This is needed for grid_area. Attribute names must have same length
shape$names <- "layer"
} else if (class(shape)[1] == "data.frame" |
class(shape)[1] == "matrix") {
## If coordinate names are found, take those columns,
## otherwise take the first 2
if (length(colnames(shape))) {
accep_cols_x <- c("L*N", "X")
accep_cols_y <- c("L*T", "Y", "BREITE")
sum_col_match <- sum(sapply(c(accep_cols_x, accep_cols_y), grepl,
toupper(colnames(shape)) ))
if (sum_col_match >= 2) {
x_col_match <- which(sapply(
lapply(accep_cols_x, grepl, toupper(colnames(shape))),
any))
y_col_match <- which(sapply(
lapply(accep_cols_y, grepl, toupper(colnames(shape))),
any))
x_col_index <- which(grepl(accep_cols_x[x_col_match],
toupper(colnames(shape))))
y_col_index <- which(grepl(accep_cols_y[y_col_match],
toupper(colnames(shape))))
pltm <- shape[, c(x_col_index[1], y_col_index[1])]
} else {
col_numeric <- which(sapply(shape[1, ], is.numeric))
pltm <- shape[, col_numeric]
}
} else {
col_numeric <- which(sapply(shape[1, ], is.numeric))
pltm <- shape[, col_numeric]
}
colnames(pltm) <- c("x", "y")
pltm <- data.frame(pltm)
shape <- st_cast(st_combine(
sf::st_as_sf(pltm, coords = c("x","y"))$geometry), "POLYGON")
if (!missing(proj)) st_crs(shape) <- 3035
}
return(shape)
}
#' @title Transform Winddata
#' @name windata_format
#' @description Helper Function, which transforms winddata to an acceptable
#' format
#'
#' @export
#'
#' @param df The wind data with speeds, direction and optionally a probability
#' column. If not assigned, it will be calculated
#'
#' @family Helper Functions
#' @return A list of windspeed and probabilities
#'
#' @examples \donttest{
#' wind_df <- data.frame(ws = c(12, 30, 45),
#' wd = c(0, 90, 150),
#' probab = 30:32)
#' windata_format(wind_df)
#'
#' wind_df <- data.frame(speed = c(12, 30, 45),
#' direction = c(90, 90, 150)
#' ,probab = c(10, 20, 60)
#' )
#' windata_format(wind_df)
#'
#' wind_df <- data.frame(speed = c(12, 30, 45),
#' direction = c(400, 90, 150)
#' )
#' windata_format(wind_df)
#' }
windata_format <- function(df) {
wind_df <- data.frame(df)
if (!all(colnames(wind_df) %in% c("ws", "wd"))) {
# Assume that we've been given a wind_df frame.
# Lets find the correct columns
if (length(colnames(wind_df)) &&
all(!colnames(wind_df) %in% c("X1", "X2", "X3")) ) {
accep_speed <- c("SPEED", "GESCH", "V", "WS")
accep_direc <- c("DIR", "RICHT", "WD")
accep_proba <- c("PRO", "WAHR")
sum_col_match <- sum(sapply(c(accep_speed, accep_direc, accep_proba),
grepl, toupper(colnames(wind_df)) ))
if (sum_col_match >= 2) {
speed_match <- which(sapply(
lapply(accep_speed, grepl, toupper(colnames(wind_df))),
any))
direc_match <- which(sapply(
lapply(accep_direc, grepl, toupper(colnames(wind_df))),
any))
probab_match <- which(sapply(
lapply(accep_proba, grepl, toupper(colnames(wind_df))),
any))
speed_index <- which(grepl(accep_speed[speed_match],
toupper(colnames(wind_df))))
direc_index <- which(grepl(accep_direc[direc_match],
toupper(colnames(wind_df))))
if (length(probab_match) != 0) {
probab_index <- which(grepl(accep_proba[probab_match],
toupper(colnames(wind_df))))
wind_df[, c(speed_index[1], direc_index[1], probab_index[1])]
colnames(wind_df) <- c("ws", "wd", "probab")
} else {
wind_df[, c(speed_index[1], direc_index[1])]
colnames(wind_df) <- c("ws", "wd")
}
} else {
col_numeric <- which(sapply(wind_df[1, ], is.numeric))
wind_df <- wind_df[, col_numeric]
colnames(wind_df) <- c("ws", "wd")
}
} else {
col_numeric <- which(sapply(wind_df[1, ], is.numeric))
wind_df <- wind_df[, col_numeric]
if (length(colnames(wind_df)) == 2) {
colnames(wind_df) <- c("ws", "wd")
} else {
colnames(wind_df) <- c("ws", "wd", "probab")
}
}
}
wind_df$wd <- round(wind_df$wd, 0)
wind_df$wd <- round(wind_df$wd / 100, 1) * 100
## If no probabilites are given, assign uniform distributed ones.
if (anyNA(colnames(wind_df))) {
which(is.na(colnames(wind_df)))
colnames(wind_df)[3] <- "probab"
}
if (any(names(wind_df) == "probab") == FALSE) {
wind_df$probab <- 100 / nrow(wind_df)
}
## Checks if all the sum of possibility is 100
if (sum(wind_df$probab) != 100) {
wind_df$probab <- wind_df$probab * (100 / sum(wind_df$probab))
}
## Checks if duplicated wind directions are at hand
if (any(duplicated(wind_df$wd))) {
for (i in 1:length(wind_df[duplicated(wind_df$wd) == FALSE, 1]) ) {
## Get duplicated direction rows
temp <- wind_df[wind_df$wd == wind_df[duplicated(
wind_df$wd) == FALSE, ][i, "wd"], ]
## Sum up speed and probability
temp$ws <- sum(temp$ws * (temp$probab / sum(temp$probab)))
temp$probab <- sum(temp$probab * (temp$probab / sum(temp$probab)))
## Assign new/uniwue windspeed and probablity per direction
wind_df[wind_df$wd == wind_df[duplicated(
wind_df$wd) == FALSE, ][i, "wd"], ]$ws <- round(temp$ws, 2)[1]
wind_df[wind_df$wd == wind_df[duplicated(
wind_df$wd) == FALSE, ][i, "wd"], ]$probab <- round(temp$probab, 2)[1]
}
}
## Delete duplicated direction rows
wind_df <- wind_df[!duplicated(wind_df$wd) == TRUE, ]
## Order by direction
wind_df <- wind_df[with(wind_df, order(wd)), ]
## Sum up probabilites to 100% again
if (sum(wind_df$probab) != 100) {
wind_df$probab <- wind_df$probab * (100 / sum(wind_df$probab))
}
probabDir <- wind_df$probab
if (any(wind_df$wd > 360)) {
wind_df[wind_df$wd > 360, "wd"] <- wind_df[wind_df$wd > 360, "wd"] - 360
}
wind_df <- as.matrix(wind_df)
winddata <- list(wind_df, probabDir)
return(winddata)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.