| Title: | Data Set and Helper Functions for Wind Farm Layout Optimization Problems |
|---|---|
| Description: | Provides a convenient data set, a set of helper functions, and a benchmark function for economically (profit) driven wind farm layout optimization. This enables researchers in the field of the NP-hard (non-deterministic polynomial-time hard) problem of wind farm layout optimization to focus on their optimization methodology contribution and also provides a realistic benchmark setting for comparability among contributions. See Croonenbroeck, Carsten & Hennecke, David (2020) <doi:10.1016/j.energy.2020.119244>. |
| Authors: | Carsten Croonenbroeck [aut, cre], David Hennecke [ctb] |
| Maintainer: | Carsten Croonenbroeck <[email protected]> |
| License: | GPL-3 |
| Version: | 1.9 |
| Built: | 2024-12-14 06:53:05 UTC |
| Source: | CRAN |
Convenience function that downloads the larger data set covering the entire area of Germany and saves it to the specified directory. This data set can replace the smaller sub sample data set that is part of this package.
AcquireData(Folder)AcquireData(Folder)
Folder |
must be a character string containing the folder location to where the file will be saved. For instance, use |
This function will interactively lead the user through the downloading process and also gives advice on how to use the larger data set. Make sure that R can write to the specified directory.
AcquireData returns nothing.
Carsten Croonenbroeck
Use FarmVars for the settings object and FarmData for the smaller, built-in data set that is part of this package.
## Not run: AcquireData(tempdir()) # Will download the data file to the specified directory. ## End(Not run)## Not run: AcquireData(tempdir()) # Will download the data file to the specified directory. ## End(Not run)
Implements a set of computations and constants to compute air density as a function of altitude (a.s.l.), temperature and relative humidity.
AirDensity(Altitude = 0, Temperature = 20, phi = 0.76)AirDensity(Altitude = 0, Temperature = 20, phi = 0.76)
Altitude |
altitude in meters above sea level. |
Temperature |
temperature in degrees Celsius. |
phi |
relative humidity, must be within [0, 1]. |
The function first computes air pressure at user provided target altitude based on the barometric formula. Then, saturation vapor pressure is computed using the Magnus formula. With both, the density of air is computed afterward.
AirDensity returns a single value representing air density in kg per cubic meter.
This function returns valid values for negative altitudes, but only for temperature values between -45 and +60 degrees Celsius.
Carsten Croonenbroeck
WindspeedLog and WindspeedHellmann to compute wind speed at target heights using two slightly different established models.
AirDensity() AirDensity(300, 25, 0.6) AirDensity(0, 0, 0) # Standard conditions according to DIN 1343:1990, # returns the expected value of 1.292.AirDensity() AirDensity(300, 25, 0.6) AirDensity(0, 0, 0) # Standard conditions according to DIN 1343:1990, # returns the expected value of 1.292.
The partial Jensen wake model requires the overlap area of the rotor disc at the downwind location with the wake disc. This function computes it.
Area(a, b, d)Area(a, b, d)
a |
must be a single value. Provide the radius of the first circle in meters. |
b |
must be a single value. Provide the radius of the second circle in meters. |
d |
must be a single value. Provide the distance between the two circles' centers. |
If a turbine is downwind another turbine, the wake cone of that upwind turbine may only partially cover the rotor disc of the second turbine. For the partial wake model it is necessary to compute the covered area.
Area returns covered area in square meters.
Carsten Croonenbroeck
JensenTrapezoid to check whether there are wake effects present. FarmVars for the data object.
Area(60, 40, 50) # Returns 2930.279.Area(60, 40, 50) # Returns 2930.279.
A function that returns the yearly installation costs for a given set of turbines (provide x and y for the turbines' locations). In its present form it only returns the 'UnitCost' value from the FarmVars settings object per turbine.
Cost(x, y)Cost(x, y)
x |
can be a single value or a numeric vector of values, contains the 'x' location(s) of turbines. |
y |
can be a single value or a numeric vector of values, contains the 'y' location(s) of turbines. |
Note that x and y should both be of length n, i.e. the numbers of values they contain should match the number of turbines in the current wind farm layout problem.
This function is a stub and can and should be replaced by something reasonable in an actual wind farm layout problem.
Cost returns a vector of values. The number of elements matches the length of x and y.
Carsten Croonenbroeck
Profit to see where to use Cost, Yield for a similar function for yearly yield.
## Returns a vector of two, c(100000, 100000). Cost(c(0.5, 0.7), c(0.2, 0.3)) ## Replace the function by another function ## also called 'Cost', embedded in environment e. ## Also, see the vignette. ## Not run: e$Cost <- function(x, y) #x, y \in R^n { retVal <- rep(e$FarmVars$UnitCost, min(length(x), length(y))) retVal[x > 0.5] <- retVal[x > 0.5] * 2 return(retVal) } set.seed(1357) NumTurbines <- 4 # For example. Result <- pso::psoptim(par = runif(NumTurbines * 2), fn = Profit, lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result rm(Cost, envir = e) ## End(Not run)## Returns a vector of two, c(100000, 100000). Cost(c(0.5, 0.7), c(0.2, 0.3)) ## Replace the function by another function ## also called 'Cost', embedded in environment e. ## Also, see the vignette. ## Not run: e$Cost <- function(x, y) #x, y \in R^n { retVal <- rep(e$FarmVars$UnitCost, min(length(x), length(y))) retVal[x > 0.5] <- retVal[x > 0.5] * 2 return(retVal) } set.seed(1357) NumTurbines <- 4 # For example. Result <- pso::psoptim(par = runif(NumTurbines * 2), fn = Profit, lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result rm(Cost, envir = e) ## End(Not run)
Environment that contains variables object FarmVars and, if downloaded, full data set FarmData.
ee
e contains data and variables, see FarmVars and FarmData.
e returns the environment for FarmVars and FarmData.
Carsten Croonenbroeck
AcquireData for downloading the full data set.
e ## Gives the environment object.e ## Gives the environment object.
This list object contains four matrices covering adjusted yield, wind speed, wind direction and standard deviations of wind directions in Germany.
FarmDataFarmData
A list object containing seven matrices, each containing 25 x 25 values at a raster resolution of 200 x 200 m (note that for the larger data set downloadable using AcquireData, each matrix contains 4400 x 3250 values):$AdjustedYield: Yield is average annual energy production (AEP). According to FGW technical guidelines, AEP is adjusted due to different location qualities to obtain a better guess at the marketable energy output. Interpret these values directly as 'megawatt hours per year'.$WindSpeed: Average wind speeds in meters per second.$WindDirection: Average wind directions in degrees (azimuth system).$SDDirection: Standard deviations of wind directions in degrees (azimuth system).$Elevation: Terrain relief (elevation in meters).$Slope: Terrain slope in degrees.$SlopeDirection: Direction of hillside in degrees.
If the full mode data set is present, it is loaded into the environment e. The built-in data set FarmData contains matrices of dimension 25 x 25, while the full data set e$FarmData consists of 4400 x 3250 matrices.
Carsten Croonenbroeck
David Hennecke
DWD Climate Data Center (CDC): ftp://opendata.dwd.de/climate_environment/CDC/grids_germany/multi_annual/wind_parameters/resol_200x200/
FGW Technical Guidelines: https://wind-fgw.de/shop/technical-guidelines/?lang=en
Croonenbroeck, C. & Hennecke, D. Does the German renewable energy act provide a fair incentive system for onshore wind power? - A simulation analysis, Energy Policy, 2020, 114, 111663, 1-15, https://doi.org/10.1016/j.enpol.2020.111663
# 'Profit' uses this data set internally: NumTurbines <- 4 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result)# 'Profit' uses this data set internally: NumTurbines <- 4 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result)
This list object collects all necessary variables for wind farm layout optimizations within this package. The object is loaded upon loading the package, embedded into a user writeable environment e and provides a set of default values.
FarmVarsFarmVars
A list object containing:$UnitCost: This contains yearly costs per turbine in actual currency. Mostly, this will be installation costs, maintenance cost and others. If, for example, total turbine installation costs are 2 mio. EUR and the projected life span of the turbine is 20 years, then the yearly cost is roughly 100,000 EUR, which is also the default value for this item. Refine the number at will to encompass annuity (taking interest effects into consideration), salvage, or the mentioned maintenance cost.$Price: This denotes the average sale price of electricity produced at the wind farm per megawatt hour. Defaults to 100, since 100 EUR per megawatt hour is a good (still rough) estimate in the European electricity market.$StartPoint: Denotes the point at which to start to 'cut out' a sample field from the data set. Since the full data set covers the entire region of Germany at a 200 x 200 meters resolution, the entire area is too big for a wind farm. For an actual wind farm layout optimization problem, a smaller region is indicated. This package will conveniently provide a square area of 5 km x 5 km in size. Note: Concerning the larger data set covering the entire area of Germany (conveniently downloadable using AcquireData), this area con be found at point 2000:2000 (full data covering a range of 4400 x 3250 points) and from there, spans a five kilometers square, i.e. covers the raster at 2000:2024, 2000:2024. Change the default value of 1 to 2000 to find that area or to something else for cross checking your algorithm with data covering a different area. However, for comparability, evaluate your algorithm result at the default of 2000 or using the built-in data set.$Width: Denotes the width (and, since the test area is a square, also the height) of the test area. Defaults to 25, which leads to a 5 km x 5 km square as a sample wind field, as the raster resolution is 200 x 200 m.$MeterMinDist: Denotes the minimum distance from one turbine to another in meters. Defaults to 500 meters, which is a rather conservative value.$EndPoint: The complement to '$StartPoint'. Is computed based on '$StartPoint' and '$Width' as follows: FarmVars$StartPoint + FarmVars$Width - 1. By default, $StartPoint = 1 and $Width = 25, so $EndPoint = 1 + 25 - 1 = 25.$MeterWidth: Contains the width (and, since the test area is a square, also the height) of the test area in meters. Since by default, FarmVars$Width = 25 and the raster resolution is 200 x 200 m, this defaults to $MeterWidth = 200 * 25 = 5000 meters (5 km).$MinDist: Contains the minimum distance in problem space. Since the problem space in a convenient unit square, this defaults to $MinDist = MeterMinDist / MeterWidth = 500 / 5000 = 0.1.$z: Contains the hub height of the turbine type under investigation. Default value is 100 meters.$z0: Contains the terrain's roughness length required for Jensen's wake model. In most studies, a default value of 0.1 is agreed to be a good guess for onshore wind farms, so this is the default value here.$r0: Contains the rotor radius for the turbine type under investigation. For instance, a Vestas V90 turbine has a rotor diameter of 90 meters (hence the name), so the radius is 45 mesters, which is also the default value here.Partial:
Selects whether to incorporate partial Jensen wake or not. Defaults to TRUE.$BenchmarkSolution: Contains the best setup for the standard benchmark field at a task of placing 20 turbines. This result has been generated by optimizer genoud() from package rgenoud. It should serve as a benchmark or starting point for improved solutions. e$FarmVars$BenchmarkSolution is taken from Croonenbroeck & Hennecke (2020). Note, however, that the computed profit is a little less than the 12,053,607 EUR reported there, since, starting from version 1.6, wflo by default takes partial Jensen wake effects into account. This does however in no way constrain the representativeness of e$FarmVars$BenchmarkSolution.
Carsten Croonenbroeck
Croonenbroeck, C. & Hennecke, D. A Comparison of Optimizers in a Unified Standard for Optimization on Wind Farm Layout Optimization, Energy, 2020, 119244, 1-15, https://doi.org/10.1016/j.energy.2020.119244
# Inspect the benchmark solution by Result <- list(par = e$FarmVars$BenchmarkSolution) Profit(Result$par) PlotResult(Result) ShowWakePenalizers(Result)# Inspect the benchmark solution by Result <- list(par = e$FarmVars$BenchmarkSolution) Profit(Result$par) PlotResult(Result) ShowWakePenalizers(Result)
Dependent on the resolution at which GenerateGauss was set to use, the resolution coordinates (array indices) cannot immediately be interpreted as being measured in meters. This function returns wind speeds from the Gauss object accoring to x, y, z simply provided in meters.
GaussWS(Gauss, x, y, z)GaussWS(Gauss, x, y, z)
Gauss |
must be an object returned by |
x |
must be a single value. Provide desired distance in meters. |
y |
must be a single value. Provide desired distance in meters. May be negative. |
z |
must be a single value. Provide desired distance in meters. |
The Gaussian wake model is loosely based on the initial contribution by Bastankhah & Porte-Agel (2014).
GaussWS returns a single number which can be considered a wind speed in the wake of a turbine at location x, y, and z.
Carsten Croonenbroeck
Bastankhah, M., & Porte-Agel, F. (2014). A new analytical model for wind-turbine wakes. Renewable Energy, 70, 116-123.
Use GenerateGauss to compute the three-dimensional tensor array object containing the wind speed data. See QuickGauss3D for the same algorithm, immediately returning the wind speed at one single point only.
## Not run: GaussWS(Gauss, 100, 1, 100) GaussWS(Gauss, 80, -40, 90) GaussWS(Gauss, 200, 40, 150) ## End(Not run)## Not run: GaussWS(Gauss, 100, 1, 100) GaussWS(Gauss, 80, -40, 90) GaussWS(Gauss, 200, 40, 150) ## End(Not run)
This function pre-computes Gaussian wind speeds and stores them in a 3D array, similar to voxels. Using that 'table', wind speeds can be looked up very quickly, which makes Gaussian wake feasible during WFLO runs.
GenerateGauss(u = 8, refHeight = 10, maxX = 500, resY = 100, resZ = 100, Verbose = TRUE)GenerateGauss(u = 8, refHeight = 10, maxX = 500, resY = 100, resZ = 100, Verbose = TRUE)
u |
measured wind speed at reference height. Will mostly be measured in meters per second. |
refHeight |
reference height in meters. This is the height at which the incoming wind speed u is measured. |
maxX |
the number of steps down the x axis for which to compute the model. |
resY |
the number of steps along the y axis for which to compute the model. Note that as y may take negative values, the resolution space should be chosen not too small, here. If, e.g., resY = 100, this means that y may take values from -50 to 50, which may be too low a resolution in some cases. |
resZ |
the number of steps up the z axis for which to compute the model. |
Verbose |
selectes whether the function displays status reports during computation, as it may take some time, dependent on the resolution setting. |
Users may choose to compute a rather fine resolution run over night and then save the returned object so it can be loaded in future sessions. The Gaussian wake model is loosely based on the initial contribution by Bastankhah & Porte-Agel (2014).
GenerateGauss returns the three-dimensional array containing wind speeds.
Note that the model assumes that along the x axis, x = 0 is the turbine location. x expands along the wind direction downwind. y denoted whether a point is 'left' or 'right' the x axis. Thus, the x-z plane is the plane along the x axis and perpendicular to the ground. The z axis is hight, starting at 0 = ground level.
Carsten Croonenbroeck
Bastankhah, M., & Porte-Agel, F. (2014). A new analytical model for wind-turbine wakes. Renewable Energy, 70, 116-123.
Use GaussWS for a convenience function to look-up the values from the returned array. See QuickGauss3D for the same algorithm, immediately returning the wind speed at one single point only.
## Not run: Gauss <- GenerateGauss(maxX = 500, resY = 1000, resZ = 1000) ## End(Not run)## Not run: Gauss <- GenerateGauss(maxX = 500, resY = 1000, resZ = 1000) ## End(Not run)
Use this function to convert degrees from the arithmetic system (0° being east, ascending counterclockwise) into the azimuth system (0° being north, ascending clockwise) and vice versa.
Geo2Ari(g)Geo2Ari(g)
g |
contains a single value degree (usually between 0° and 360°, decimal fractions allowed). |
g may contain degrees from both systems, the function turns the data into the respective other system.
Geo2Ari returns a single value of degrees.
Carsten Croonenbroeck
GetDirInfo for further degree information.
Geo2Ari(0) ## In an arithmetic system, 0° means 'east', while 'east' in ## azimuth notation is 90°. This call returns 90. Geo2Ari(90) ## In an azimuth system, 90° means 'east', while 'east' in ## arithmetic notation is 0°. This call returns 0. Geo2Ari(Geo2Ari(123)) ## Returns 123.Geo2Ari(0) ## In an arithmetic system, 0° means 'east', while 'east' in ## azimuth notation is 90°. This call returns 90. Geo2Ari(90) ## In an azimuth system, 90° means 'east', while 'east' in ## arithmetic notation is 0°. This call returns 0. Geo2Ari(Geo2Ari(123)) ## Returns 123.
As seen from a point in the wind farm, computes the angle to another point in that farm.
GetAngle(x1, y1, x2, y2)GetAngle(x1, y1, x2, y2)
x1 |
must be a single value. Provide the |
y1 |
must be a single value. Provide the |
x2 |
must be a single value. Provide the |
y2 |
must be a single value. Provide the |
From point 2's point of view, computes the angle to point 1.
GetAngle returns a single number between 0 and 360. If both points are identical, the return value is 0.
Note that this function returns an angle in arithmetic degrees notation. To convert to azimuth notation, use Geo2Ari.
Carsten Croonenbroeck
Use JensenAngle to compute the wake cone and with it, use JensenTrapezoid to see if another turbine B is in turbine A's wake.
GetAngle(0.2, 0.2, 0.1, 0.1) ## Looking from point (0.1, 0.1) at point (0.2, 0.2), the angle is 45° (arithmetic).GetAngle(0.2, 0.2, 0.1, 0.1) ## Looking from point (0.1, 0.1) at point (0.2, 0.2), the angle is 45° (arithmetic).
PlotResult.Given a point, an angle information (in arithmetic degrees), and a length information, computes start and end points of an arrow usable via the arrows function.
GetArrow(BaseX, BaseY, Degrees, Frac = 25)GetArrow(BaseX, BaseY, Degrees, Frac = 25)
BaseX |
must be a single value containing the x value of a point which later will be the center of the arrow. |
BaseY |
must be a single value containing the y value of a point which later will be the center of the arrow. |
Degrees |
must be a single value containing the desired rotation degree of the arrow. |
Frac |
must be a single value containing the length of the arrow. Default is 25 and for convenience, this parameter should in most cases be identical to FarmVars$Width. |
This function will be used internally by PlotResult.
GetArrow returns a vector of four values representing x and y for the start point and x and y for the end point of an arrow (in that order).
Carsten Croonenbroeck
Use PlotResult to visualize the optimization result. See FarmVars for the data object.
GetArrow(0.5, 0.5, 45) #At c(0.5, 0.5), generates an arrow pointing in north-eastern direction.GetArrow(0.5, 0.5, 45) #At c(0.5, 0.5), generates an arrow pointing in north-eastern direction.
For a turbine's location represented by x and y, looks up the wind direction from the matrix Dirs and the corresponding standard deviation from matrix SDs. Internally transforms coordinates of x and y from problem space (usually unit square) to the matrix space of Dirs and SDs, respectively.
GetDirInfo(x, y, Dirs, SDs)GetDirInfo(x, y, Dirs, SDs)
x |
must be a single value containing the 'x' location of a turbine in problem space. |
y |
must be a single value containing the 'y' location of a turbine in problem space. |
Dirs |
a matrix containing average yearly wind directions. Usually, the third element of the list object |
SDs |
a matrix containing average yearly wind direction standard deviations. Usually, the fourth element of the list object |
GetDirInfo returns a vector of two elements, the first being the average wind direction in degrees and the second being the corresponding standard deviation. Note that degrees are meant in the arithmetic degrees system (0° being east, ascending counterclockwise). To transform into an azimuth system (0° being north, ascending clockwise), use function Geo2Ari. Also note that wind directions are meant to denote 'where the wind is going to' rather than 'where the wind is coming from'.
Carsten Croonenbroeck
Profit to see where to use GetDirInfo, Yield for a similar function for adjusted yield. See FarmVars for the data object.
Dirs <- FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] GetDirInfo(0.5, 0.7, Dirs, SDs) ## Returns wind direction and standard deviation for the given location ## and provided the matrices given.Dirs <- FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] GetDirInfo(0.5, 0.7, Dirs, SDs) ## Returns wind direction and standard deviation for the given location ## and provided the matrices given.
This function accepts a pair longitude (decimal system, east) and latitude (decimal system, north) coordinates as well a a desired edge length in m and if valid, returns a list object similar to FarmData, but at the specified location and with the specified size.
GetFarmFromLonLat(Top, Left, EdgeLen, DoPlot = FALSE)GetFarmFromLonLat(Top, Left, EdgeLen, DoPlot = FALSE)
Top |
longitude value (decimal system, east). |
Left |
latitude value (decimal system, north). |
EdgeLen |
edge length in meters. Note that the returned farm data object is always a square area and thus, contains square matrices. |
DoPlot |
optionally plots the annual energy production (AEP) 'landscape' in the returned object using a |
Requires that the full FarmData dataset is loaded. See AcquireData on how to obtain it.
Returns a list object following the structure of FarmData, but only containing the farm area starting from the top-left point specified and as many 'tiles' required to meet the desired edge length at 200 m tiles resolution.
Carsten Croonenbroeck
See Index2GK to convert index coordinates to Gauss-Kruger coordinates, and GK2LonLat to convert Gauss-Kruger coordinates to longitude/latitude coordinates.
# This will return a farm at the specified location, edge length 5,000 m (5 km). # Requires full data set to be loaded. ## Not run: MyFarm <- GetFarmFromLonLat(51.49594, 11.58818, 5000, DoPlot = FALSE) ## End(Not run)# This will return a farm at the specified location, edge length 5,000 m (5 km). # Requires full data set to be loaded. ## Not run: MyFarm <- GetFarmFromLonLat(51.49594, 11.58818, 5000, DoPlot = FALSE) ## End(Not run)
This function accepts Gauss-Kruger (GK) coordinates that are within the FarmData space and converts them to a pair of indices that point to locations within the FarmData matrices.
GK2Index(X, Y)GK2Index(X, Y)
X |
a numeric value for the 'right' value in the GK coordinate system. |
Y |
a numeric value for the 'top' value in the GK coordinate system. |
Assumes that the GK zone is EPSG 31467. References to indices valid if the full data set is loaded, see FarmData and AcquireData.
Returns index values between 1 and 3,250 for X and between 1 and 4,400 for Y.
Carsten Croonenbroeck
Index2GK for the inverse.
GK2Index(3280000, 5230000) # Will return c(1, 4400), the lower left point. GK2Index(3929800, 6109800) # Will return c(3250, 1), the top right point.GK2Index(3280000, 5230000) # Will return c(1, 4400), the lower left point. GK2Index(3929800, 6109800) # Will return c(3250, 1), the top right point.
This function accepts a pair of Gauss-Kruger (GK) coordinates and converts them to longitude (decimal system, east) and latitude (decimal system, north) coordinates by performing re-projection.
GK2LonLat(X, Y)GK2LonLat(X, Y)
X |
a numeric value for the 'right' value in the GK coordinate system. |
Y |
a numeric value for the 'top' value in the GK coordinate system. |
Assumes that the GK zone is EPSG 31467.
Returns a vector of two values where longitude is first, latitude second.
Carsten Croonenbroeck
LonLat2GK for the inverse.
LonLat <- GK2LonLat(3702793, 5998319) # Will return c(12.09750, 54.07547).LonLat <- GK2LonLat(3702793, 5998319) # Will return c(12.09750, 54.07547).
For a turbine's location represented by x and y, looks up the elevation from the matrix Elev. Internally transforms coordinates of x and y from problem space (usually unit square) to the matrix space of Elev.
Height(x, y, Elev)Height(x, y, Elev)
x |
must be a single value containing the 'x' location of a turbine in problem space. |
y |
must be a single value containing the 'y' location of a turbine in problem space. |
Elev |
a matrix containing heights. Usually, the fifth element of the list object |
Height is a convenience function for looking up heights as required e.g. for a partial Jensen wake model, independent from the actual size of the area under investigation.
Height returns a single value, the elevation in meters.
Carsten Croonenbroeck
Profit to see where to use Height. See FarmData for the data set.
## Returns adjusted yield for the given location. P <- c(0.5868695, 0.9722714) Height(P[1], P[2], FarmData[[5]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint])## Returns adjusted yield for the given location. P <- c(0.5868695, 0.9722714) Height(P[1], P[2], FarmData[[5]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint])
PlotResult.Draws a set of arrows over an existing plot. Will be used internally by PlotResult to visualize the vector field of wind directions over a plot of the wind farm.
ImposeVectorField(xNum = 5, yNum = 5, ColMain = "black", ColBand = "white", Frac = 25, DoSDs = TRUE)ImposeVectorField(xNum = 5, yNum = 5, ColMain = "black", ColBand = "white", Frac = 25, DoSDs = TRUE)
xNum |
must be a single value containing the desired number of arrows horizontally. |
yNum |
must be a single value containing the desired number of arrows vertically. |
ColMain |
must be a single value containing the desired color for the main wind direction arrow. |
ColBand |
must be a single value containing the desired color for the secondary arrows representing the standard deviation arrows. Only used if |
Frac |
must be a single value containing the desired length information for the arrows. Will be passed to |
DoSDs |
must be |
This function will be used internally by PlotResult.ImposeVectorField requires FarmData to be present and an existing plot to impose the vector field over.
ImposeVectorField returns nothing.
Carsten Croonenbroeck
Use PlotResult to visualize the optimization result.
plot(c(0, 1), c(0, 1)) ImposeVectorField(ColBand = "red")plot(c(0, 1), c(0, 1)) ImposeVectorField(ColBand = "red")
This function accepts indices that are within the FarmData space and converts them to a pair of Gauss-Kruger (GK) coordinates.
Index2GK(X, Y)Index2GK(X, Y)
X |
a numeric value for the 'X' value in the FarmData space. |
Y |
a numeric value for the 'Y' value in the FarmData space. |
Assumes that the GK zone is EPSG 31467. References to indices are valid if the full data set is loaded, see FarmData and AcquireData.
For index values between 1 and 3,250 for X and between 1 and 4,400 for Y, returns Gauss-Kruger coordinates valid in the EPSG 31467 coordinate system.
Carsten Croonenbroeck
GK2Index for the inverse.
Index2GK(1, 4400) # Will return c(3280000, 5230000), the lower left point. Index2GK(3250, 1) # Will return c(3929800, 6109800), the top right point.Index2GK(1, 4400) # Will return c(3280000, 5230000), the lower left point. Index2GK(3250, 1) # Will return c(3929800, 6109800), the top right point.
x, computes the wake cone generated by a turbine.In Jensen's wake model, a wake cone is generated by a turbine based on the turbine's hub height z, its rotor radius r_0, the terrain's roughness length z_0 and the distance x between the turbine under investigation and a certain second point. z, r_0 and z_0 are taken from the FarmVars settings object, while x is to be provided.
JensenAngle(x)JensenAngle(x)
x |
must be a single value. Provide distance in meters. |
If a second turbine (B) is in a first turbine (A)'s wake, turbine B's power output must be penalized due to turbulence. This can be performed via Jensen's wake model, but first, it should be checked whether B is actually in A's wake. This function computes the angle of the wake cone generated by turbine A. This can be used together with function JensenTrapezoid PointInPolygon and to see if B is in the wake of A.
JensenAngle returns the (onesided) angle of the wake cone generated by a turbine for distance x.
Carsten Croonenbroeck
Jensen, N. O. (1983). A note on wind generator interaction. Roskilde: Risø National Laboratory. Risø-M, No. 2411
JensenTrapezoid to check whether there are wake effects present. FarmVars for the data object.
JensenAngle(500) ## For a distance of 500 m, the function returns a wake cone angle of 9.2233 degrees ## (given standard values for z, z_0 and r_0 in the FarmVars settings object).JensenAngle(500) ## For a distance of 500 m, the function returns a wake cone angle of 9.2233 degrees ## (given standard values for z, z_0 and r_0 in the FarmVars settings object).
x, computes the penalty factor for a turbine's wake.In Jensen's wake model, a wake cone is generated by a turbine based on the turbine's hub height z, its rotor radius r_0, the terrain's roughness length z_0 and the distance x between the turbine under investigation and a certain second point. z, r_0 and z_0 are taken from the FarmVars settings object, while x is to be provided.
JensenFactor(x)JensenFactor(x)
x |
must be a single value. Provide distance in meters. |
If a second turbine (B) is in a first turbine (A)'s wake, turbine B's power output must be penalized due to turbulence. This function computes the penalty factor for that wake.
JensenFactor returns the a single number between 0 and 1, which can immediately be used a a penalty factor.
Note that for Jensen's multiple wake model, wake penalties from several potentially influencing upwind turbines must be computed. As the factors are between 0 and 1, they can conveniently multiplied by each other afterward. If there are four turbines (A, B, C, D) inside a wind farm, you are investigating turbine D and it turns out that D is in the wake of B and C, but not A, then compute penalty for B (assume 0.7) and C (assume 0.8). Now, the penalty factors for A, B and C are 1, 0.7 and 0.8, so 1 * 0.7 * 0.8 = 0.56 is the penalty for turbine D.
Carsten Croonenbroeck
Jensen, N. O. (1983). A note on wind generator interaction. Roskilde: Risø National Laboratory. Risø-M, No. 2411
Use JensenAngle to compute the wake cone and with it, use JensenTrapezoid to see if another turbine B is in turbine A's wake. Apply JensenFactor only if this is the case.
JensenFactor(500) ## Provided that turbine B is in turbine A's wake and that the two turbines are 500 meters apart, ## turbine B must be penalized by a factor of 0.7952.JensenFactor(500) ## Provided that turbine B is in turbine A's wake and that the two turbines are 500 meters apart, ## turbine B must be penalized by a factor of 0.7952.
Provided a wind direction, a point and a downwind length, computes the corner points of a Jensen trapezoid (in the literature oftentimes called a 'cone').
JensenTrapezoid(WindDir, Point, x, Margin = TRUE)JensenTrapezoid(WindDir, Point, x, Margin = TRUE)
WindDir |
unique value, a wind direction in degrees. |
Point |
a vector of two containing x and y coordinates of a point. Usually, both will be between 0 and 1. |
x |
downwind distance of the cone in meters. |
Margin |
specified whether to add a small margin to the distance parameter (i.e. internally, compute x * 1.1). Defaults to TRUE. |
JensenTrapezoid returns a matrix of four columns (the four corner points) and two rows (x and y coordinates).
Carsten Croonenbroeck
Jensen, N. O. (1983). A note on wind generator interaction. Roskilde: Riso National Laboratory. Riso-M, No. 2411
JensenAngle to compute 'cone' angle information, PointInPolygon for a test whether a point is inside a polygon.
MyTrapezoid <- JensenTrapezoid(45, c(0.5, 0.5), 500)MyTrapezoid <- JensenTrapezoid(45, c(0.5, 0.5), 500)
This function accepts a pair longitude (decimal system, east) and latitude (decimal system, north) coordinates and converts them to Gauss-Kruger (GK) coordinates by performing re-projection.
LonLat2GK(Lon, Lat)LonLat2GK(Lon, Lat)
Lon |
longitude value (decimal system, east). |
Lat |
latitude value (decimal system, north). |
The resulting values are in GK zone EPSG 31467.
Returns a vector of two values where X is first, Y second.
Carsten Croonenbroeck
GK2LonLat for the inverse.
LonLat2GK(12.09750, 54.07548) # Will return c(3702793, 5998320).LonLat2GK(12.09750, 54.07548) # Will return c(3702793, 5998320).
This function returns the dimensionless cost based on Mosetti et al. (1994).
MosettiTurbineCost(N)MosettiTurbineCost(N)
N |
number of turbines for which cost is to be computed. |
Wind farm cost does not increase linearly with the number of turbines. Economies of scale lead to disproportionate cost increases. This function returns the dimensionless cost based on Mosetti et al. (1994).
MosettiTurbineCost returns a dimensionless cost factor.
Carsten Croonenbroeck
G. Mosetti, C. Poloni, B. Diviacco (1994). Optimization of wind turbine positioning in large windfarms by means of a genetic algorithm, Journal of Wind Engineering and Industrial Aerodynamics, 51(1), 105-116.
Cost for a plug-in cost function for profit.
MosettiTurbineCost(1) MosettiTurbineCost(20)MosettiTurbineCost(1) MosettiTurbineCost(20)
As seen from a turbine in the wind farm, computes the wake penalty factor for another turbine in that farm.
PairPenalty(x1, y1, x2, y2, Dirs, SDs)PairPenalty(x1, y1, x2, y2, Dirs, SDs)
x1 |
must be a single value. Provide the |
y1 |
must be a single value. Provide the |
x2 |
must be a single value. Provide the |
y2 |
must be a single value. Provide the |
Dirs |
a matrix containing average yearly wind directions. Usually, the third element of the list object |
SDs |
a matrix containing average yearly wind direction standard deviations. Usually, the fourth element of the list object |
First, this function uses GetAngle to compute the angle between the two points provided, as seen from point 2's point of view. It then obtains the wind direction at point 2 using GetDirInfo. After that, the distance between the two points is computed. With it, the wake cone is computed using JensenAngle to check whether point 2 is in point 1's wake using JensenTrapezoid. If that is the case, JensenFactor is used to compute the penalty factor.
Note that the penalty is the deduction to wind speed. It applies to wind power by its third power, so the user is responsible to take it to its cube himself. Profit does that automatically internally.
PairPenalty returns a single number between 0 and 1. If point 2 is not in the wake of point 1, the function returns 1.
Carsten Croonenbroeck
Use JensenFactor to see how this function operates. See FarmVars for the data object.
Dirs <- FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] PairPenalty(0.9, 0.8, 0.6, 0.9, Dirs, SDs) ## Weak wake penalty PairPenalty(0.1, 0.1, 0.6, 0.9, Dirs, SDs) ## No wake penaltyDirs <- FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] PairPenalty(0.9, 0.8, 0.6, 0.9, Dirs, SDs) ## Weak wake penalty PairPenalty(0.1, 0.1, 0.6, 0.9, Dirs, SDs) ## No wake penalty
If a rotor disc of a first turbine is only partially covered by a wake cone of a second turbine, the penalty must be adjusted accordingly. This function returns partial penalties, if partial wake applies, or full penalty, if not.
PartialJensen(x2, y2, x1, y1, Dirs, SDs, Elev, Z = NULL, DrawTop = FALSE, DrawFront = FALSE, DrawSide = FALSE)PartialJensen(x2, y2, x1, y1, Dirs, SDs, Elev, Z = NULL, DrawTop = FALSE, DrawFront = FALSE, DrawSide = FALSE)
x1 |
must be a single value. Provide the |
y1 |
must be a single value. Provide the |
x2 |
must be a single value. Provide the |
y2 |
must be a single value. Provide the |
Dirs |
a matrix containing average yearly wind directions. Usually, the third element of the list object |
SDs |
a matrix containing average yearly wind direction standard deviations. Usually, the fourth element of the list object |
Elev |
a matrix containing elevations. Usually, the fifth element of the list object |
Z |
accepts a vector of two representing the heights of the turbines at the two points in meters. If |
DrawTop |
If |
DrawFront |
If |
DrawSide |
If |
The function first checks whether there is partial coverage. If so, it adjusts the penalty internally given by PairPenalty and if not, returns the full penalty. Therefore, the function is imposed over an existing Jensen model and refines it.
Note that the penalty is the deduction to wind speed. It applies to wind power by its third power, so the user is responsible to take it to its cube himself. Profit does that automatically internally.
PartialJensen returns a penality between 0 and 1.
Carsten Croonenbroeck
Frandsen, S. (1992). On the wind speed reduction in the center of large clusters of wind turbines. Journal of Wind Engineering and Industrial Aerodynamics, 39(1-3), pp. 251-265, https://doi.org/10.1016/0167-6105(92)90551-K.
JensenTrapezoid to check whether there are wake effects present. FarmVars for the data object. PairPenalty for the non-partial wake penalty.
P1 <- c(0.5868695, 0.9722714) P2 <- c(0.4827957, 0.9552658) if (exists("FarmData", envir = e, inherits = FALSE)) { Dirs <- e$FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- e$FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] Elev <- e$FarmData[[5]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] } else { Dirs <- FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] Elev <- FarmData[[5]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] } ## First, compute the non-partial penalty: Penal <- PairPenalty(P2[1], P2[2], P1[1], P1[2], Dirs, SDs) ## Then, correct it for partial coverage: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev) ## Now draw a top view: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev, DrawTop = TRUE) ## Now draw a front view: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev, DrawFront = TRUE) ## Note that the downwind point is somewhat elevated, seems 'right' ## from the upwind point of view, and far away (rotor disc seems smaller). ## Now draw a side view: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev, DrawSide = TRUE) ## For elevation, see Height(P1[1], P1[2], Elev) # (upwind point) ## and Height(P2[1], P2[2], Elev) # (downwind point) ## Instead, for the next example, the downwind point is closer ## to the upwind point (larger impression of rotor disc), 'left' ## of it and lower in terms of elevation: x1 <- 0.5 y1 <- 0.5 x2 <- 0.45 y2 <- 0.478 Penal <- PairPenalty(x2, y2, x1, y1, Dirs, SDs) Penal2 <- PartialJensen(x2, y2, x1, y1, Dirs, SDs, Elev, DrawTop = TRUE) Penal2 <- PartialJensen(x2, y2, x1, y1, Dirs, SDs, Elev, DrawFront = TRUE) ## For elevation, see Height(x1, y1, Elev) # (upwind point) ## and Height(x2, y2, Elev) # (downwind point)P1 <- c(0.5868695, 0.9722714) P2 <- c(0.4827957, 0.9552658) if (exists("FarmData", envir = e, inherits = FALSE)) { Dirs <- e$FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- e$FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] Elev <- e$FarmData[[5]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] } else { Dirs <- FarmData[[3]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] SDs <- FarmData[[4]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] Elev <- FarmData[[5]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] } ## First, compute the non-partial penalty: Penal <- PairPenalty(P2[1], P2[2], P1[1], P1[2], Dirs, SDs) ## Then, correct it for partial coverage: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev) ## Now draw a top view: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev, DrawTop = TRUE) ## Now draw a front view: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev, DrawFront = TRUE) ## Note that the downwind point is somewhat elevated, seems 'right' ## from the upwind point of view, and far away (rotor disc seems smaller). ## Now draw a side view: Penal2 <- PartialJensen(P2[1], P2[2], P1[1], P1[2], Dirs, SDs, Elev, DrawSide = TRUE) ## For elevation, see Height(P1[1], P1[2], Elev) # (upwind point) ## and Height(P2[1], P2[2], Elev) # (downwind point) ## Instead, for the next example, the downwind point is closer ## to the upwind point (larger impression of rotor disc), 'left' ## of it and lower in terms of elevation: x1 <- 0.5 y1 <- 0.5 x2 <- 0.45 y2 <- 0.478 Penal <- PairPenalty(x2, y2, x1, y1, Dirs, SDs) Penal2 <- PartialJensen(x2, y2, x1, y1, Dirs, SDs, Elev, DrawTop = TRUE) Penal2 <- PartialJensen(x2, y2, x1, y1, Dirs, SDs, Elev, DrawFront = TRUE) ## For elevation, see Height(x1, y1, Elev) # (upwind point) ## and Height(x2, y2, Elev) # (downwind point)
This function draws the adjusted yields of the wind farm under investigation using image, superimposes a contour plot using contour and arrows for the wind directions using ImposeVectorField and then draws the points for the turbines' locations (aligned to their respective raster grid centers).
PlotResult(Result, ImageData = NULL, DoLabels = FALSE, Labels = "IDs")PlotResult(Result, ImageData = NULL, DoLabels = FALSE, Labels = "IDs")
Result |
must be an optimization result as returned by an optimizer such as |
ImageData |
a matrix containing the data for the background, processed via |
DoLabels |
a boolean that indicates whether labels should be plotted next to the points. Defaults to FALSE. |
Labels |
a vector of length n = N / 2, can be numeric values or strings. Defines the labels to be shown if |
For maximum convenience and compatibility with numeric optimizers of most kinds, this function expects nothing but the usual optimization result list. This, however, requires that the FarmData dataset as well as the FarmVars object is present defining additional settings.
PlotResult returns nothing.
Carsten Croonenbroeck
Use Profit to obtain an optimization result.
#Will not provide a very good result NumTurbines <- 4 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result)#Will not provide a very good result NumTurbines <- 4 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result)
Uses the famous Jensen test to check whether a provided point is inside a polygon, the latter defined by a set of points.
PointInPolygon(PointMat, TestPoint)PointInPolygon(PointMat, TestPoint)
PointMat |
a 2 x 4 matrix, four corner points of each x and y coordinates. |
TestPoint |
a vector of two containing x and y coordinates of a point to test. |
PointInPolygon returns TRUE if the point is inside the polygon, or FALSE, else.
Carsten Croonenbroeck
Jensen, N. O. (1983). A note on wind generator interaction. Roskilde: Risø National Laboratory. Risø-M, No. 2411
JensenTrapezoid to compute a Jensen trapezoid.
set.seed(1357) Angle <- runif(n = 1, min = 0, max = 360) Point <- runif(n = 2, min = 0, max = 1) Dist <- rnorm(n = 1, mean = 200, sd = 200) MyTrapezoid <- JensenTrapezoid(Angle, Point, Dist) plot(x = MyTrapezoid[1, ], y = MyTrapezoid[2, ], xlab = "", ylab = "") lines(x = c(MyTrapezoid[1, 1], MyTrapezoid[1, 2]), y = c(MyTrapezoid[2, 1], MyTrapezoid[2, 2])) lines(x = c(MyTrapezoid[1, 2], MyTrapezoid[1, 3]), y = c(MyTrapezoid[2, 2], MyTrapezoid[2, 3])) lines(x = c(MyTrapezoid[1, 3], MyTrapezoid[1, 4]), y = c(MyTrapezoid[2, 3], MyTrapezoid[2, 4])) lines(x = c(MyTrapezoid[1, 4], MyTrapezoid[1, 1]), y = c(MyTrapezoid[2, 4], MyTrapezoid[2, 1])) NumTest <- 50 xTest <- runif(n = NumTest, min = min(MyTrapezoid[1, ]), max = max(MyTrapezoid[1, ])) yTest <- runif(n = NumTest, min = min(MyTrapezoid[2, ]), max = max(MyTrapezoid[2, ])) for (i in 1:NumTest) { ThisPoint <- c(xTest[i], yTest[i]) if (PointInPolygon(MyTrapezoid, ThisPoint)) { points(ThisPoint[1], ThisPoint[2], pch = 16, col = "green") } else { points(ThisPoint[1], ThisPoint[2], pch = 16, col = "red") } }set.seed(1357) Angle <- runif(n = 1, min = 0, max = 360) Point <- runif(n = 2, min = 0, max = 1) Dist <- rnorm(n = 1, mean = 200, sd = 200) MyTrapezoid <- JensenTrapezoid(Angle, Point, Dist) plot(x = MyTrapezoid[1, ], y = MyTrapezoid[2, ], xlab = "", ylab = "") lines(x = c(MyTrapezoid[1, 1], MyTrapezoid[1, 2]), y = c(MyTrapezoid[2, 1], MyTrapezoid[2, 2])) lines(x = c(MyTrapezoid[1, 2], MyTrapezoid[1, 3]), y = c(MyTrapezoid[2, 2], MyTrapezoid[2, 3])) lines(x = c(MyTrapezoid[1, 3], MyTrapezoid[1, 4]), y = c(MyTrapezoid[2, 3], MyTrapezoid[2, 4])) lines(x = c(MyTrapezoid[1, 4], MyTrapezoid[1, 1]), y = c(MyTrapezoid[2, 4], MyTrapezoid[2, 1])) NumTest <- 50 xTest <- runif(n = NumTest, min = min(MyTrapezoid[1, ]), max = max(MyTrapezoid[1, ])) yTest <- runif(n = NumTest, min = min(MyTrapezoid[2, ]), max = max(MyTrapezoid[2, ])) for (i in 1:NumTest) { ThisPoint <- c(xTest[i], yTest[i]) if (PointInPolygon(MyTrapezoid, ThisPoint)) { points(ThisPoint[1], ThisPoint[2], pch = 16, col = "green") } else { points(ThisPoint[1], ThisPoint[2], pch = 16, col = "red") } }
This is the bread-and-butter function to this package. It takes a set of points (turbine locations) and based on the adjusted yields in the field specified via the FarmVars settings object, checks whether the layout is valid, computes the Jensen penalties, generates the total marketable power production of this layout, takes cost and sale price into consideration and finally computes the farm's economic profit. Note, however, that since most numeric optimizers by default operate as a minimizer, Profit returns the negative profit, i.e. a negative number for positive profit values and a positive number if cost is greater than revenue.
Profit(X)Profit(X)
X |
must be a numeric vector containing an even number of values, at least two. If the length of the vector is N, then n = N / 2 is the number of points. The first values 1,...,n are interpreted as x coordinates and the subsequent n + 1,...,N values are the y coordinates. |
For maximum convenience and compatibility with numeric optimizers of most kinds, this function expects nothing but the problem vector (the set of points) as a parameter. This, however, requires that the FarmData dataset as well as the FarmVars object is present defining additional settings.
As Profit returns the negative profit for compatibility with minimizers, reverse the sign for actual profit values.Profit requires the optimizer to accept box constraints, as Profit can not compute values outside the wind farm boundary (the function's domain). If your optimizer does not accept box constraints, embed Profit into a wrapper function that returns the sum of costs if it least one point is outside the boundary.
Profit returns a single number. The result is the negative profit for any valid setting of points and consequently, the sum of all costs for invalid settings.
Carsten Croonenbroeck
Calls Yield and Cost internally. Use PlotResult to visualize the optimization result.
#Will not provide a very good result NumTurbines <- 4 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result) ########################################## #Will provide a somewhat better result #Necessary to install pso ## Not run: NumTurbines <- 4 Result <- pso::psoptim(par = runif(NumTurbines * 2), fn = Profit, lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result) ## End(Not run) ########################################## #Simple wrapper function for optimizers not accepting box constraints NumTurbines <- 4 lower <- rep(0, NumTurbines * 2) upper <- rep(1, NumTurbines * 2) Wrapper <- function(X) { xSel <- seq(from = 1, to = length(X) - 1, by = 2) x <- X[xSel] y <- X[xSel + 1] if (any(x < lower) | any(x > upper) | any(y < lower) | any(y > upper)) { return(sum(rep(e$FarmVars$UnitCost, length(x)))) } return(Profit(X)) } ## Not run: set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Wrapper, method = "SANN") Result PlotResult(Result) ## End(Not run)#Will not provide a very good result NumTurbines <- 4 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result) ########################################## #Will provide a somewhat better result #Necessary to install pso ## Not run: NumTurbines <- 4 Result <- pso::psoptim(par = runif(NumTurbines * 2), fn = Profit, lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result PlotResult(Result) ## End(Not run) ########################################## #Simple wrapper function for optimizers not accepting box constraints NumTurbines <- 4 lower <- rep(0, NumTurbines * 2) upper <- rep(1, NumTurbines * 2) Wrapper <- function(X) { xSel <- seq(from = 1, to = length(X) - 1, by = 2) x <- X[xSel] y <- X[xSel + 1] if (any(x < lower) | any(x > upper) | any(y < lower) | any(y > upper)) { return(sum(rep(e$FarmVars$UnitCost, length(x)))) } return(Profit(X)) } ## Not run: set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Wrapper, method = "SANN") Result PlotResult(Result) ## End(Not run)
This function takes an optimizer result (or simple vector object) and computes the profit contribution for each point.
ProfitContributors(Result)ProfitContributors(Result)
Result |
must be an optimization result as returned by an optimizer such as |
Neglecting a possibly invalid setup (caused by at least one point), computes the profit contributions. If the setup provided is valid, the sum of contributions is identical to the (absulute value of the) returned value of Profit.
ProfitContributors returns a matrix of two columns and n rows, whith n being the number of turbines. The first column is a sequence of 1:n, representing the turbine IDs, while the second column contains the actual profit contributions.
Carsten Croonenbroeck
See ValidSetup to see how a setup is categorized as valid or not. Use PlotResult to visualize the optimization result.
#Prints a result and uses the profit contributions as labels. NumTurbines <- 4 set.seed(1235) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) MyLabels <- ProfitContributors(Result) MyLabels #PlotResult(Result, DoLabels = TRUE, Labels = MyLabels[, 2]) # Given a valid setup, this should be TRUE. #identical(abs(Profit(Result$par)), sum(MyLabels[, 2]))#Prints a result and uses the profit contributions as labels. NumTurbines <- 4 set.seed(1235) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) MyLabels <- ProfitContributors(Result) MyLabels #PlotResult(Result, DoLabels = TRUE, Labels = MyLabels[, 2]) # Given a valid setup, this should be TRUE. #identical(abs(Profit(Result$par)), sum(MyLabels[, 2]))
This function computes Gaussian wind speeds.
QuickGauss3D(x, y, z, u = 8, refHeight = 10)QuickGauss3D(x, y, z, u = 8, refHeight = 10)
x |
must be a single value. Provide distance in meters. |
y |
must be a single value. Provide distance in meters. |
z |
must be a single value. Provide distance in meters. |
u |
must be a single value. Provide wind speed in meters per second. |
refHeight |
must be a single value. Provide reference height (the height at which wind speed u was measured) in meters. |
The Gaussian wake model is loosely based on the initial contribution by Bastankhah & Porte-Agel (2014).
QuickGauss3D returns a single number which can be considered a wind speed in the wake of a turbine at location x, y, and z.
Note that the model assumes that along the x axis, x = 0 is the turbine location. x expands along the wind direction downwind. y denoted whether a point is 'left' or 'right' the x axis. Thus, the x-z plane is the plane along the x axis and perpendicular to the ground. The z axis is hight, starting at 0 = ground level.
Carsten Croonenbroeck
Bastankhah, M., & Porte-Agel, F. (2014). A new analytical model for wind-turbine wakes. Renewable Energy, 70, 116-123.
Use GenerateGauss to compute the three-dimensional tensor array object containing the wind speed data. Use GaussWS for a convenience function to look-up the values from the returned array.
QuickGauss3D(100, 1, 100) QuickGauss3D(200, -40, 120) QuickGauss3D(50, 40, 70)QuickGauss3D(100, 1, 100) QuickGauss3D(200, -40, 120) QuickGauss3D(50, 40, 70)
After optimization, this function draws a reduced 'field' and shows which points cause wake penalties on which points. 'Causers' are displayed in grey, 'sufferers', i.e. points that are in (possibly multiple) wakes of (several) causers are displayed in gold.
ShowWakePenalizers(Result, Cones = TRUE, All = FALSE, VectorField = TRUE, ImposeMST = FALSE, NS = FALSE, NSOnly = FALSE, Exaggerate = TRUE, DoubleLine = TRUE, Frames = 100, Alpha = 0.5, MaxContrast = TRUE, Soften = TRUE)ShowWakePenalizers(Result, Cones = TRUE, All = FALSE, VectorField = TRUE, ImposeMST = FALSE, NS = FALSE, NSOnly = FALSE, Exaggerate = TRUE, DoubleLine = TRUE, Frames = 100, Alpha = 0.5, MaxContrast = TRUE, Soften = TRUE)
Result |
must be an optimization result as returned by an optimizer such as |
Cones |
triggers whether the Jensen cones are to be displayed as well. Defaults to |
All |
triggers whether all Jensen cones are to be displayed, instead of just causers and sufferers. Defaults to |
VectorField |
controls whether or not to impose the wind directions vector field (arrows). Defaults to |
ImposeMST |
if |
NS |
if |
NSOnly |
if |
Exaggerate |
if |
DoubleLine |
applies only if |
Frames |
number of frames to pre-compute for the Navier-Stokes simulation. 100 is a good value. Less makes the computation faster, but more inaccurate. More than 100 in only necessary for very low wind speeds. Thus, defaults to |
Alpha |
The Navier-Stokes wake simulation is imposed semi-transparent over the usual image. This controls the transparency. Alpha = 0 means transparent, Alpha = 1 means opaque. Defaults to |
MaxContrast |
increases the contrast to maximum (this is equivalent to histogram stretching). Defaults to |
Soften |
applies a Gaussian soften filter to the Navier-Stokes image to have a more homogeneous visualization. Defaults to |
Enables the researcher to inspect the optimization result in an a-posteriori analysis of Jensen's wake model. Parts of the Navier-Stokes code are inspired by Stam (2003).
ShowWakePenalizers invisibly returns a square matrix (dimension: 'number of turbines') of wake penalty pairs: Columns are sufferers, rows are causers. For each row/column pair unequal to one, the turbine in row i sheds wake on the turbine in column j. The actual numbers are those returned by JensenFactor.
Carsten Croonenbroeck
Stam, S. (2003). Real-Time Fluid Dynamics for Games. https://www.researchgate.net/publication/2560062_Real-Time_Fluid_Dynamics_for_Games, 1-17
Use PlotResult to visualize the optimization result.
## Not run: #Will show that turbine 1 sheds wake on turbine 5. NumTurbines <- 8 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) ShowWakePenalizers(Result) ############### #This may take some time. A progress bar is displayed. Result = list(par = e$FarmVars$BenchmarkSolution) ShowWakePenalizers(Result, All = TRUE, NS = TRUE, VectorField = FALSE, Alpha = 0.7) ############### #Generates a few frames of an animation. Result <- list(par = e$FarmVars$BenchmarkSolution) MakeFrames <- function(Frames) { for (i in 100:(100 + Frames - 1)) { bmp(file = paste(i, ".bmp", sep = "")) ShowWakePenalizers(Result, All = TRUE, NS = TRUE, Frames = i) dev.off() } } system.time(MakeFrames(30)) #May take an hour. ## End(Not run)## Not run: #Will show that turbine 1 sheds wake on turbine 5. NumTurbines <- 8 set.seed(1357) Result <- optim(par = runif(NumTurbines * 2), fn = Profit, method = "L-BFGS-B", lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) ShowWakePenalizers(Result) ############### #This may take some time. A progress bar is displayed. Result = list(par = e$FarmVars$BenchmarkSolution) ShowWakePenalizers(Result, All = TRUE, NS = TRUE, VectorField = FALSE, Alpha = 0.7) ############### #Generates a few frames of an animation. Result <- list(par = e$FarmVars$BenchmarkSolution) MakeFrames <- function(Frames) { for (i in 100:(100 + Frames - 1)) { bmp(file = paste(i, ".bmp", sep = "")) ShowWakePenalizers(Result, All = TRUE, NS = TRUE, Frames = i) dev.off() } } system.time(MakeFrames(30)) #May take an hour. ## End(Not run)
This function helps with analzing an already computed turbine setup. As Profit (together with an optimizer) finds an optimal setup given actual wind speeds in that area, the question arises whether that setup is 'robust' against changing wind directions. If for a given setup, there is no or only little wake dependence (analyze e.g. with ShowWakePenalizers), this may change if wind direction changes. Given a different wind direction, a situation may come up at which there is an upwind turbine and a local set of, say, two, three, or even more downwind turbines. As these downwind turbines are in the wake of the upwind turbines ('receiving' only reduced wind speeds), it may be beneficial to shut down the upwind turbine in order to have the downwind turbine receive full wind, thus possibly generating greater profit.SwitchProfile expects a 'bit' pattern of turbines (so-called 'profiles') at which '1' means that a turbine is to be used, while '0' means, the turbine is to be shut down. For a total of four tubines in a setup, profile c(1, 1, 1, 1) means that all turbines are used, while profile c(1, 0, 0, 1) means that turbines 1 and 4 are used, turbines 2 and 3 are to be shut down. If a wind farm is 'robust' against wind direction changes, it should return profiles of all ones for all wind directions.
In the examples we show a short code sample that uses SwitchProfile to loop through 360 integer wind direction degrees and optimizes actual profiles. It is up to the researcher to use alternative optimizers in order to find (possibly) better profiles.
SwitchProfile(On, Result)SwitchProfile(On, Result)
On |
must be a numeric vector of length |
Result |
the actual setup as an optimizer result. Must be a list that contains the turbine locations is a slot "$par". For example, use |
For compatibility with non-integer solution optimizers, SwitchProfile expects not only values equal to 0 or 1, but will assume that values below 0.5 are meant to be 0, and values >= 0.5 are coerced to 1.
SwitchProfile returns a single number. The result is the negative profit for the profile given. The result can be optimized in order to find the maximum profit profile, which should be identical to a profile of all ones if the setup is robust against wind direction changes.
Carsten Croonenbroeck
Use Profit for the target function used internally.
Result <- list(par = e$FarmVars$BenchmarkSolution) Profile1 <- rep(1, times = 20) #Will be identical to Profit(Result$par): SwitchProfile(Profile1, Result) Profile2 <- Profile1 Profile2[8] <- 0 #Disable turbine 8. #Returns farm profit if turbine 8 is off: SwitchProfile(Profile2, Result) ########################################## #Warning, this will overwrite e$FarmData. If not dispensable, #backup before use. #Computation may be slow. ## Not run: Result <- list(par = e$FarmVars$BenchmarkSolution) e$FarmData <- FarmData n <- length(e$FarmVars$BenchmarkSolution) / 2 Profiles <- matrix(ncol = n, nrow = 360) Dom <- cbind(rep(0, n), rep(1, n)) for (j in 1:360) { e$FarmData$WindDirection <- matrix(data = j, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) message(paste("Processing wind direction: ", j, " degrees.", sep = "")) flush.console() Result <- rgenoud::genoud(fn = SwitchProfile, nvars = n, starting.values = runif(n), Domains = Dom, boundary.enforcement = 2, MemoryMatrix = FALSE, print.level = 0, Result = Result) On <- Result$par On[On < 0.5] <- 0 On[On >= 0.5] <- 1 Profiles[j, ] <- On } #Now 'Profiles' is a matrix of 360 rows (degrees) and n columns. #If the setup is robust against wind direction changes, all values #in the columns should be ones. If there are zeros left, double #check these profiles (compare to full turbines setup) manually, #e.g. using for (i in 1:nrow(Profiles)) { if (any(Profiles[i, ] == 0)) { e$FarmData$WindDirection <- matrix(data = i, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) print(paste("Is unrestricted profit (", abs(Profit(e$FarmVars$BenchmarkSolution)), ") still greater than 'profile profit' (", abs(SwitchProfile(Profiles[i, ], Result)), ")? ", abs(Profit(e$FarmVars$BenchmarkSolution)) >= abs(SwitchProfile(Profiles[i, ], Result)), ".", sep = "")) } } #Is absolute profile profit actually greater than the absolute value #function Profit() returns? If yes, then it is beneficial to turn off #turbines (use profiles) if wind comes from the given angle. ## End(Not run) ########################################## #As an alternative to using an optimizer to find the optimal profile, #all possible profiles can be 'looped through' deterministically using #the following few lines of code. #Warning, this will overwrite e$FarmData. If not dispensable, #backup before use. #Be advised, computation may be very slow, even if run on modern #machines and even though using parallelization. ## Not run: Result <- list(par = e$FarmVars$BenchmarkSolution) e$FarmData <- FarmData ComputeProfile <- function(i) { e$FarmData$WindDirection <- matrix(data = i, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) BestProfit <- n * e$FarmVars$UnitCost BestProfile <- rep(0, n) for (j in 0:MaxNum) { ThisProfile <- as.numeric(rev(intToBits(j)))[(32 - (n - 1)):32] if (length(which(ThisProfile == 1)) > 1) { ThisProfit <- SwitchProfile(ThisProfile, Result) } else { ThisProfit <- n * e$FarmVars$UnitCost } if (ThisProfit < BestProfit) { BestProfit <- ThisProfit BestProfile <- ThisProfile } } return(BestProfile) } n <- length(e$FarmVars$BenchmarkSolution) / 2 MaxNum <- (2 ^ n) - 1 NumCores <- parallel::detectCores() - 1 cl <- snow::makeCluster(NumCores) ParallelProfile <- new.env() ParallelProfile$FarmVars <- e$FarmVars ParallelProfile$FarmData <- e$FarmData snow::clusterExport(cl, list = ls()) snow::clusterExport(cl, list = ls(envir = as.environment("package:wflo"))) doSNOW::registerDoSNOW(cl) iterations <- 360 pb <- txtProgressBar(max = iterations, style = 3) progress <- function(n) setTxtProgressBar(pb, n) opts <- list(progress = progress) `%dopar%` <- foreach::`%dopar%` AllProfiles <- foreach::foreach(i = 1:iterations, .options.snow = opts) %dopar% ComputeProfile(i) close(pb) snow::stopCluster(cl) ## End(Not run) #Afterward, check the resulting list of profiles for remaining zeros as above. ########################################## #Thinking the thought of wind direction dependent layouts further, #Profit() can be used to optimize a farm setup that takes all wind #direction wake patterns into account during the optimization #stage already. #Warning, this will overwrite \code{e$FarmData}. If not dispensable, #backup before use. #Computation may be slow. ## Not run: e$FarmData <- FarmData SumProfit <- function(X) { MySum <- as.numeric() for (i in 1:360) { e$FarmData$WindDirection <- matrix(data = i, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) MySum[i] <- Profit(X) } MySum <- sum(MySum) return(MySum) } n <- 20 Dom <- cbind(rep(0, 2 * n), rep(1, 2 * n)) Result <- rgenoud::genoud(fn = SumProfit, nvars = 2 * n, starting.values = runif(2 * n), Domains = Dom, boundary.enforcement = 2, MemoryMatrix = FALSE, print.level = 1) ## End(Not run)Result <- list(par = e$FarmVars$BenchmarkSolution) Profile1 <- rep(1, times = 20) #Will be identical to Profit(Result$par): SwitchProfile(Profile1, Result) Profile2 <- Profile1 Profile2[8] <- 0 #Disable turbine 8. #Returns farm profit if turbine 8 is off: SwitchProfile(Profile2, Result) ########################################## #Warning, this will overwrite e$FarmData. If not dispensable, #backup before use. #Computation may be slow. ## Not run: Result <- list(par = e$FarmVars$BenchmarkSolution) e$FarmData <- FarmData n <- length(e$FarmVars$BenchmarkSolution) / 2 Profiles <- matrix(ncol = n, nrow = 360) Dom <- cbind(rep(0, n), rep(1, n)) for (j in 1:360) { e$FarmData$WindDirection <- matrix(data = j, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) message(paste("Processing wind direction: ", j, " degrees.", sep = "")) flush.console() Result <- rgenoud::genoud(fn = SwitchProfile, nvars = n, starting.values = runif(n), Domains = Dom, boundary.enforcement = 2, MemoryMatrix = FALSE, print.level = 0, Result = Result) On <- Result$par On[On < 0.5] <- 0 On[On >= 0.5] <- 1 Profiles[j, ] <- On } #Now 'Profiles' is a matrix of 360 rows (degrees) and n columns. #If the setup is robust against wind direction changes, all values #in the columns should be ones. If there are zeros left, double #check these profiles (compare to full turbines setup) manually, #e.g. using for (i in 1:nrow(Profiles)) { if (any(Profiles[i, ] == 0)) { e$FarmData$WindDirection <- matrix(data = i, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) print(paste("Is unrestricted profit (", abs(Profit(e$FarmVars$BenchmarkSolution)), ") still greater than 'profile profit' (", abs(SwitchProfile(Profiles[i, ], Result)), ")? ", abs(Profit(e$FarmVars$BenchmarkSolution)) >= abs(SwitchProfile(Profiles[i, ], Result)), ".", sep = "")) } } #Is absolute profile profit actually greater than the absolute value #function Profit() returns? If yes, then it is beneficial to turn off #turbines (use profiles) if wind comes from the given angle. ## End(Not run) ########################################## #As an alternative to using an optimizer to find the optimal profile, #all possible profiles can be 'looped through' deterministically using #the following few lines of code. #Warning, this will overwrite e$FarmData. If not dispensable, #backup before use. #Be advised, computation may be very slow, even if run on modern #machines and even though using parallelization. ## Not run: Result <- list(par = e$FarmVars$BenchmarkSolution) e$FarmData <- FarmData ComputeProfile <- function(i) { e$FarmData$WindDirection <- matrix(data = i, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) BestProfit <- n * e$FarmVars$UnitCost BestProfile <- rep(0, n) for (j in 0:MaxNum) { ThisProfile <- as.numeric(rev(intToBits(j)))[(32 - (n - 1)):32] if (length(which(ThisProfile == 1)) > 1) { ThisProfit <- SwitchProfile(ThisProfile, Result) } else { ThisProfit <- n * e$FarmVars$UnitCost } if (ThisProfit < BestProfit) { BestProfit <- ThisProfit BestProfile <- ThisProfile } } return(BestProfile) } n <- length(e$FarmVars$BenchmarkSolution) / 2 MaxNum <- (2 ^ n) - 1 NumCores <- parallel::detectCores() - 1 cl <- snow::makeCluster(NumCores) ParallelProfile <- new.env() ParallelProfile$FarmVars <- e$FarmVars ParallelProfile$FarmData <- e$FarmData snow::clusterExport(cl, list = ls()) snow::clusterExport(cl, list = ls(envir = as.environment("package:wflo"))) doSNOW::registerDoSNOW(cl) iterations <- 360 pb <- txtProgressBar(max = iterations, style = 3) progress <- function(n) setTxtProgressBar(pb, n) opts <- list(progress = progress) `%dopar%` <- foreach::`%dopar%` AllProfiles <- foreach::foreach(i = 1:iterations, .options.snow = opts) %dopar% ComputeProfile(i) close(pb) snow::stopCluster(cl) ## End(Not run) #Afterward, check the resulting list of profiles for remaining zeros as above. ########################################## #Thinking the thought of wind direction dependent layouts further, #Profit() can be used to optimize a farm setup that takes all wind #direction wake patterns into account during the optimization #stage already. #Warning, this will overwrite \code{e$FarmData}. If not dispensable, #backup before use. #Computation may be slow. ## Not run: e$FarmData <- FarmData SumProfit <- function(X) { MySum <- as.numeric() for (i in 1:360) { e$FarmData$WindDirection <- matrix(data = i, ncol = e$FarmVars$Width, nrow = e$FarmVars$Width) MySum[i] <- Profit(X) } MySum <- sum(MySum) return(MySum) } n <- 20 Dom <- cbind(rep(0, 2 * n), rep(1, 2 * n)) Result <- rgenoud::genoud(fn = SumProfit, nvars = 2 * n, starting.values = runif(2 * n), Domains = Dom, boundary.enforcement = 2, MemoryMatrix = FALSE, print.level = 1) ## End(Not run)
For a set of turbine locations represented by x and y, checks whether all possible pairs satisfy the minimum distance criterion.
ValidSetup(x, y)ValidSetup(x, y)
x |
must be a numeric vector of at least two values, contains the 'x' location(s) of turbines. |
y |
must be a numeric vector of at least two values, contains the 'y' location(s) of turbines. |
ValidSetup returns TRUE if all pairs of turbines are at least as far away from each other as 'MinDist' from the FarmVars settings object requests, or FALSE, else.
Carsten Croonenbroeck
ValidSetup(c(0.5, 0.7), c(0.7, 0.9)) ## Returns TRUE if FarmVars$MinDist is at its standard value (0.1).ValidSetup(c(0.5, 0.7), c(0.7, 0.9)) ## Returns TRUE if FarmVars$MinDist is at its standard value (0.1).
This package makes two contributions to the Wind Farm Layout Optimization (WFLO) research branch: First, it provides a convenient and realistic data set of high resolution and accuracy encompassing wind speeds, wind directions, standard deviations of wind directions, and (adjusted) yields for the entire area of Germany. Second, it supplies a set of helper functions and a benchmark function for economically (profit) driven wind farm layout optimization. This enables researchers in the field of the np-hard problem of wflo to focus on their optimization methodology contribution and also provides a realistic benchmark setting for comparability among contributions.
Carsten Croonenbroeck
David Hennecke
For a wind speed at given height, returns a scaled wind speed at some different height. Often used to obtain wind speed at hub height.
WindspeedHellmann(v0, HH = 100, refHeight = 10, alpha = 1 / 7)WindspeedHellmann(v0, HH = 100, refHeight = 10, alpha = 1 / 7)
v0 |
wind speed (in meters per second) at reference height. |
HH |
the height (in meters) at which wind speed is desired. |
refHeight |
reference height (in meters). The height at which the actual wind speed (v0) was measured. |
alpha |
power law parameter. Usually set to 1/7 for onshore sites. |
This function simply implements
WindspeedHellmann returns a wind speed (in meters per second) at the desired height.
Carsten Croonenbroeck
WindspeedLog for a different way to scale wind speeds to heights.
WindspeedLog(v0 = 6, HH = 80, z0 = 0.1, refHeight = 20) WindspeedHellmann(v0 = 6, HH = 80, refHeight = 20)WindspeedLog(v0 = 6, HH = 80, z0 = 0.1, refHeight = 20) WindspeedHellmann(v0 = 6, HH = 80, refHeight = 20)
For a wind speed at given height, returns a scaled wind speed at some different height. Often used to obtain wind speed at hub height.
WindspeedLog(v0, HH = 100, z0 = 0.1, refHeight = 10)WindspeedLog(v0, HH = 100, z0 = 0.1, refHeight = 10)
v0 |
wind speed (in meters per second) at reference height. |
HH |
the height (in meters) at which wind speed is desired. |
z0 |
roughness length (in meters). Usually set to 0.1 m for onshore sites. |
refHeight |
reference height (in meters). The height at which the actual wind speed (v0) was measured. |
This function simply implements
Note that this way to scale wind speeds to certain heights is frequently considered deprecated in the literature. Use WindspeedHellmann instead.
WindspeedLog returns a wind speed (in meters per second) at the desired height.
Carsten Croonenbroeck
WindspeedHellmann for a different way to scale wind speeds to heights.
WindspeedLog(v0 = 6, HH = 80, z0 = 0.1, refHeight = 20) WindspeedHellmann(v0 = 6, HH = 80, refHeight = 20)WindspeedLog(v0 = 6, HH = 80, z0 = 0.1, refHeight = 20) WindspeedHellmann(v0 = 6, HH = 80, refHeight = 20)
For a turbine's location represented by x and y, looks up the (adjusted) yield from the matrix Adj. Internally transforms coordinates of x and y from problem space (usually unit square) to the matrix space of Adj.
Yield(x, y, Adj)Yield(x, y, Adj)
x |
must be a single value containing the 'x' location of a turbine in problem space. |
y |
must be a single value containing the 'y' location of a turbine in problem space. |
Adj |
a matrix containing adjusted yields. Usually, the first element of the list object |
Adjusted yields are the projected yearly average yields dependent on wind speed, hub height and other settings at each point in the raster data. Annual Energy Production (AEP) at a specific location, weighted by a location quality correction factor, produces adjusted yields. This adjustment returns a better guess on the marketable yield at a specific point. For details on the data, see the data set description to this package.
Note that Profit internally multiplies the outcome of Yield by e$FarmVars$Price to obtain revenue. Users who replace the function by e$Yield need to provide that manually, if revenue is desired.
Yield returns a single value.
Carsten Croonenbroeck
Profit to see where to use Yield, Cost for a similar function for yearly cost. FarmData for the data set.
## Returns adjusted yield for the given location. Adj <- FarmData[[1]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] Yield(0.5, 0.7, Adj) ## Replace the function by another function ## also called 'Yield', embedded in environment e. ## Also, see the vignette. ## Not run: e$Yield <- function(x, y, AEP) #x, y \in R { return(x + y) } set.seed(1357) NumTurbines <- 4 # For example. Result <- pso::psoptim(par = runif(NumTurbines * 2), fn = Profit, lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result rm(Yield, envir = e) ## End(Not run)## Returns adjusted yield for the given location. Adj <- FarmData[[1]][e$FarmVars$StartPoint:e$FarmVars$EndPoint, e$FarmVars$StartPoint:e$FarmVars$EndPoint] Yield(0.5, 0.7, Adj) ## Replace the function by another function ## also called 'Yield', embedded in environment e. ## Also, see the vignette. ## Not run: e$Yield <- function(x, y, AEP) #x, y \in R { return(x + y) } set.seed(1357) NumTurbines <- 4 # For example. Result <- pso::psoptim(par = runif(NumTurbines * 2), fn = Profit, lower = rep(0, NumTurbines * 2), upper = rep(1, NumTurbines * 2)) Result rm(Yield, envir = e) ## End(Not run)