From e416754156820c7785a704a0471d11796520b07e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 16 Apr 2021 16:57:44 +0200 Subject: [PATCH 01/20] Added functions (and power curves) to compute wind power density and wind capacity factor --- R/WindCapacityFactor.R | 15 ++++++++++++++ R/WindPowerDensity.R | 8 ++++++++ R/zzz.R | 45 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 68 insertions(+) create mode 100644 R/WindCapacityFactor.R create mode 100644 R/WindPowerDensity.R diff --git a/R/WindCapacityFactor.R b/R/WindCapacityFactor.R new file mode 100644 index 0000000..f6dc9f3 --- /dev/null +++ b/R/WindCapacityFactor.R @@ -0,0 +1,15 @@ + +WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/III", "III")) { + IEC_class <- match.arg(IEC_class) + pc_files <- c( + "I" = "Enercon_E70_2.3MW.txt", + "I/II" = "Gamesa_G80_2.0MW.txt", + "II" = "Gamesa_G87_2.0MW.txt", + "II/III" = "Vestas_V100_2.0MW.txt", + "III" = "Vestas_V110_2.0MW.txt" + ) + pc_file <- system.file("power_curves", pc_files[IEC_class], package="CSIndicators", mustWork=T) + pc <- read_pc(pc_file) + cf <- wind2CF(wind, pc) + return(cf) +} diff --git a/R/WindPowerDensity.R b/R/WindPowerDensity.R new file mode 100644 index 0000000..70407fd --- /dev/null +++ b/R/WindPowerDensity.R @@ -0,0 +1,8 @@ + +#======================= +# Wind power density +# default: 1.225, the standard density of air at 15 ºC and 1013.25 hPa +#======================= +WindPowerDensity <- function(wind, ro = 1.225) +{ return(0.5 * ro * wind^3) +} diff --git a/R/zzz.R b/R/zzz.R index 33660aa..365d362 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -23,3 +23,48 @@ } return(position) } + + +#======================= +# Read a powercurve file +# Create the approximation function +# TODO: Document file format +#======================= +read_pc <- function(file) { + pc <- list() + + # Read pc points + pc$points <- rbind(c(0,0), read.delim(file,comment.char="#")) + + # Create an approximating function + pc$fun <- approxfun(pc$points$WindSpeed,pc$points$Power,method="linear",yleft=NA,yright=0) + + # Read other turbine characteristics (using perl's grep) + attr <- strsplit(trimws(system(paste("perl -e 'open FH,\"",file,"\";while(){@parts= /^# (.+): (.+) /;print \"@parts \";}'",sep=""),intern=T)),"\\s+") + attr <- matrix(unlist(attr), ncol = 2, byrow = T) + pc$attr <- as.list(attr[, 2]) + names(pc$attr) <- attr[, 1] + pc$attr$Filename <- file + + # TODO: check required attributes! + pc$attr$RatedPower <- as.numeric(pc$attr$RatedPower) + + return(pc) +} + +#======================= +# Evaluate the linear piecewise approximation function with the wind speed inputs to get wind power +#======================= +wind2power <- function(wind, pc) +{ power <- pc$fun(wind) + return(power) +} + +#======================= +# Convert wind to power, and divide by rated power to obtain Capacity Factor values +#======================= +wind2CF <- function(wind, pc) +{ power <- wind2power(wind, pc) + CF <- power / pc$attr$RatedPower + return(CF) +} -- GitLab From 3e878e183042544a33c3a85decbe9cdc4dde1761 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 16 Apr 2021 16:59:04 +0200 Subject: [PATCH 02/20] Adding five power curves --- NAMESPACE | 2 ++ inst/power_curves/Enercon_E70_2.3MW.txt | 38 ++++++++++++++++++++++ inst/power_curves/Gamesa_G80_2.0MW.txt | 37 ++++++++++++++++++++++ inst/power_curves/Gamesa_G87_2.0MW.txt | 36 +++++++++++++++++++++ inst/power_curves/Vestas_V100_2.0MW.txt | 42 +++++++++++++++++++++++++ inst/power_curves/Vestas_V110_2.0MW.txt | 41 ++++++++++++++++++++++++ 6 files changed, 196 insertions(+) create mode 100644 inst/power_curves/Enercon_E70_2.3MW.txt create mode 100644 inst/power_curves/Gamesa_G80_2.0MW.txt create mode 100644 inst/power_curves/Gamesa_G87_2.0MW.txt create mode 100644 inst/power_curves/Vestas_V100_2.0MW.txt create mode 100644 inst/power_curves/Vestas_V110_2.0MW.txt diff --git a/NAMESPACE b/NAMESPACE index 09b5709..08ad46b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,6 +19,8 @@ export(SelectPeriodOnDates) export(Threshold) export(TotalSpellTimeExceedingThreshold) export(TotalTimeExceedingThreshold) +export(WindPowerDensity) +export(WindCapacityFactor) import(multiApply) importFrom(ClimProjDiags,Subset) importFrom(s2dv,Reorder) diff --git a/inst/power_curves/Enercon_E70_2.3MW.txt b/inst/power_curves/Enercon_E70_2.3MW.txt new file mode 100644 index 0000000..da4be61 --- /dev/null +++ b/inst/power_curves/Enercon_E70_2.3MW.txt @@ -0,0 +1,38 @@ +#--------------------------- +# TURBINE CHARACTERISTICS +#--------------------------- +# Name: Enercon_E70_2.3MW [string] +# Manufacturer: Enercon [string] +# Diameter: 71 m +# CutIn: 2 m/s +# CutOut: 25 m/s +# ReCutIn: unknown m/s +# RatedSpeed: 16 m/s +# RatedPower: 2310 kW +# IECClass: Ia [Ia/IIa/IIIa/S] +# Control: pitch [pitch/stall/active_stall/flaps] +# HubHeights: 57,64,74,85,98,113 m +# +#--------------------------- +# POWER CURVE +# Density: 1.225 kg/m^3 +#--------------------------- +WindSpeed Power +1.0 0 +2.0 2 +3.0 18 +4.0 56 +5.0 127 +6.0 240 +7.0 400 +8.0 626 +9.0 892 +10.0 1223 +11.0 1590 +12.0 1900 +13.0 2080 +14.0 2230 +15.0 2300 +16.0 2310 +25.0 2310 +25.5 0 diff --git a/inst/power_curves/Gamesa_G80_2.0MW.txt b/inst/power_curves/Gamesa_G80_2.0MW.txt new file mode 100644 index 0000000..0b3c0bf --- /dev/null +++ b/inst/power_curves/Gamesa_G80_2.0MW.txt @@ -0,0 +1,37 @@ +#--------------------------- +# TURBINE CHARACTERISTICS +#--------------------------- +# Name: Gamesa_G80_2.0MW [string] +# Manufacturer: Gamesa [string] +# Diameter: 80 m +# CutIn: 4 m/s +# CutOut: 25 m/s +# ReCutIn: unknown m/s +# RatedSpeed: 17 m/s +# RatedPower: 2000 kW +# IECClass: Ia/IIa [Ia/IIa/IIIa/S] +# Control: pitch [pitch/stall/active_stall/flaps] +# HubHeights: 60,67,78,100 m +# +#--------------------------- +# POWER CURVE +# Density: 1.225 kg/m^3 +#--------------------------- +WindSpeed Power +3.0 0 +4.0 66 +5.0 152 +6.0 280 +7.0 457 +8.0 690 +9.0 978 +10.0 1296 +11.0 1598 +12.0 1818 +13.0 1935 +14.0 1980 +15.0 1995 +16.0 1999 +17.0 2000 +25.0 2000 +25.5 0 diff --git a/inst/power_curves/Gamesa_G87_2.0MW.txt b/inst/power_curves/Gamesa_G87_2.0MW.txt new file mode 100644 index 0000000..eefee50 --- /dev/null +++ b/inst/power_curves/Gamesa_G87_2.0MW.txt @@ -0,0 +1,36 @@ +#--------------------------- +# TURBINE CHARACTERISTICS +#--------------------------- +# Name: Gamesa_G87_2.0MW [string] +# Manufacturer: Gamesa [string] +# Diameter: 87 m +# CutIn: 4 m/s +# CutOut: 25 m/s +# ReCutIn: unknown m/s +# RatedSpeed: 16 m/s +# RatedPower: 2000 kW +# IECClass: IIa [Ia/IIa/IIIa/S] +# Control: pitch [pitch/stall/active_stall/flaps] +# HubHeights: 67,78,90,100 m +# +#--------------------------- +# POWER CURVE +# Density: 1.225 kg/m^3 +#--------------------------- +WindSpeed Power +3.0 0 +4.0 79 +5.0 181 +6.0 335 +7.0 550 +8.0 832 +9.0 1175 +10.0 1530 +11.0 1816 +12.0 1963 +13.0 1988 +14.0 1996 +15.0 1999 +16.0 2000 +25.0 2000 +25.5 0 diff --git a/inst/power_curves/Vestas_V100_2.0MW.txt b/inst/power_curves/Vestas_V100_2.0MW.txt new file mode 100644 index 0000000..462aeaa --- /dev/null +++ b/inst/power_curves/Vestas_V100_2.0MW.txt @@ -0,0 +1,42 @@ +#--------------------------- +# TURBINE CHARACTERISTICS +#--------------------------- +# Name: Vestas_V100_2.0MW [string] +# Manufacturer: Vestas [string] +# Diameter: 100 m +# CutIn: 3 m/s +# CutOut: 20 m/s +# ReCutIn: 18 m/s +# RatedSpeed: 12 m/s +# RatedPower: 2000 kW +# IECClass: IIa/IIIa [Ia/IIa/IIIa/S] +# Control: pitch [pitch/stall/active_stall/flaps] +# +#--------------------------- +# POWER CURVE +# Density: 1.225 kg/m^3 +# noise_mode: 0 +#--------------------------- +WindSpeed Power +2.5 0 +3 13 +3.5 51 +4 107 +4.5 175 +5 253 +5.5 346 +6 454 +6.5 584 +7 738 +7.5 912 +8 1109 +8.5 1321 +9 1538 +9.5 1734 +10 1873 +10.5 1951 +11 1984 +11.5 1995 +12 2000 +20 2000 +20.5 0 diff --git a/inst/power_curves/Vestas_V110_2.0MW.txt b/inst/power_curves/Vestas_V110_2.0MW.txt new file mode 100644 index 0000000..5e7657c --- /dev/null +++ b/inst/power_curves/Vestas_V110_2.0MW.txt @@ -0,0 +1,41 @@ +#--------------------------- +# TURBINE CHARACTERISTICS +#--------------------------- +# Name: Vestas_V110_2.0MW [string] +# Manufacturer: Vestas [string] +# Diameter: 110 m +# CutIn: 3 m/s +# CutOut: 20 m/s +# ReCutIn: unknown m/s +# RatedSpeed: 11.5 m/s +# RatedPower: 2000 kW +# IECClass: IIIa [Ia/IIa/IIIa/S] +# Control: pitch [pitch/stall/active_stall/flaps] +# HubHeights: 80,95,125 m +# +#--------------------------- +# POWER CURVE +# Density: 1.225 kg/m^3 +#--------------------------- +WindSpeed Power +2.5 0 +3.0 23 +3.5 80 +4.0 140 +4.5 223 +5.0 314 +5.5 422 +6.0 549 +6.5 703 +7.0 900 +7.5 1123 +8.0 1347 +8.5 1555 +9.0 1775 +9.5 1907 +10.0 1972 +10.5 1997 +11.0 1999 +11.5 2000 +20.0 2000 +20.5 0 -- GitLab From 778a5b99d6bc9a0daa1b9b6e4d8c0868a5143ccd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 16 Apr 2021 17:07:21 +0200 Subject: [PATCH 03/20] Correcting R style issues --- R/zzz.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/zzz.R b/R/zzz.R index 365d362..636943d 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -34,13 +34,13 @@ read_pc <- function(file) { pc <- list() # Read pc points - pc$points <- rbind(c(0,0), read.delim(file,comment.char="#")) + pc$points <- rbind(c(0, 0), read.delim(file, comment.char = "#")) # Create an approximating function - pc$fun <- approxfun(pc$points$WindSpeed,pc$points$Power,method="linear",yleft=NA,yright=0) + pc$fun <- approxfun(pc$points$WindSpeed, pc$points$Power, method = "linear", yleft = NA, yright = 0) # Read other turbine characteristics (using perl's grep) - attr <- strsplit(trimws(system(paste("perl -e 'open FH,\"",file,"\";while(){@parts= /^# (.+): (.+) /;print \"@parts \";}'",sep=""),intern=T)),"\\s+") + attr <- strsplit(trimws(system(paste("perl -e 'open FH,\"", file, "\";while(){@parts= /^# (.+): (.+) /;print \"@parts \";}'", sep = ""), intern = T)), "\\s+") attr <- matrix(unlist(attr), ncol = 2, byrow = T) pc$attr <- as.list(attr[, 2]) names(pc$attr) <- attr[, 1] -- GitLab From 99108f686c1cf79564c6c224beb6d4cdc63c0b12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 16 Apr 2021 17:07:56 +0200 Subject: [PATCH 04/20] Correcting R style issues --- R/WindPowerDensity.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/WindPowerDensity.R b/R/WindPowerDensity.R index 70407fd..21ed6cb 100644 --- a/R/WindPowerDensity.R +++ b/R/WindPowerDensity.R @@ -3,6 +3,6 @@ # Wind power density # default: 1.225, the standard density of air at 15 ºC and 1013.25 hPa #======================= -WindPowerDensity <- function(wind, ro = 1.225) -{ return(0.5 * ro * wind^3) +WindPowerDensity <- function(wind, ro = 1.225) { + return(0.5 * ro * wind^3) } -- GitLab From 46e49bd4a82b19fbef9708bd9fdd2290c933bc42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 16 Apr 2021 17:08:17 +0200 Subject: [PATCH 05/20] Correcting R style issues --- R/WindCapacityFactor.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/WindCapacityFactor.R b/R/WindCapacityFactor.R index f6dc9f3..15517bd 100644 --- a/R/WindCapacityFactor.R +++ b/R/WindCapacityFactor.R @@ -8,7 +8,7 @@ WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/III", "II/III" = "Vestas_V100_2.0MW.txt", "III" = "Vestas_V110_2.0MW.txt" ) - pc_file <- system.file("power_curves", pc_files[IEC_class], package="CSIndicators", mustWork=T) + pc_file <- system.file("power_curves", pc_files[IEC_class], package = "CSIndicators", mustWork = T) pc <- read_pc(pc_file) cf <- wind2CF(wind, pc) return(cf) -- GitLab From 1562f0d262d475d2a80748d9e92359166538c95c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 10:40:12 +0200 Subject: [PATCH 06/20] Read rated power from pc values, to avoid perl dependence --- R/zzz.R | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/R/zzz.R b/R/zzz.R index 636943d..35cd5a6 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -28,7 +28,6 @@ #======================= # Read a powercurve file # Create the approximation function -# TODO: Document file format #======================= read_pc <- function(file) { pc <- list() @@ -39,15 +38,8 @@ read_pc <- function(file) { # Create an approximating function pc$fun <- approxfun(pc$points$WindSpeed, pc$points$Power, method = "linear", yleft = NA, yright = 0) - # Read other turbine characteristics (using perl's grep) - attr <- strsplit(trimws(system(paste("perl -e 'open FH,\"", file, "\";while(){@parts= /^# (.+): (.+) /;print \"@parts \";}'", sep = ""), intern = T)), "\\s+") - attr <- matrix(unlist(attr), ncol = 2, byrow = T) - pc$attr <- as.list(attr[, 2]) - names(pc$attr) <- attr[, 1] - pc$attr$Filename <- file - - # TODO: check required attributes! - pc$attr$RatedPower <- as.numeric(pc$attr$RatedPower) + # Get the rated power from the power values + pc$attr$RatedPower <- max(pc$points$Power) return(pc) } -- GitLab From 086704d13391ee49f2635c3b41f6876b31199c58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 10:40:38 +0200 Subject: [PATCH 07/20] added documentation to the functions --- R/WindCapacityFactor.R | 18 +++++++++++++++++- R/WindPowerDensity.R | 20 +++++++++++++++----- 2 files changed, 32 insertions(+), 6 deletions(-) diff --git a/R/WindCapacityFactor.R b/R/WindCapacityFactor.R index 15517bd..dba82f9 100644 --- a/R/WindCapacityFactor.R +++ b/R/WindCapacityFactor.R @@ -1,4 +1,20 @@ - +#'Wind capacity factor +#' +#'@author Llorenç Lledó, \email{llledo@bsc.es} +#'@description Wind capacity factor computes the wind power generated by a specific wind turbine model under specific wind speed conditions, and expresses it as a fraction of the rated capacity (i.e. maximum power) of the turbine. +#'@description It is computed by means of a tabular power curve that relates wind speed to power output. The tabular values are interpolated with a linear piecewise approximating function to obtain a smooth power curve. Five different power curves that span different IEC classes can be selected (see below). +#'@references Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F. J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 +#'@references International Standard IEC 61400-1 (third ed.) (2005) +#' +#'@param wind a multidimensional array, vector or scalar with instantaneous wind speeds expressed in m/s. +#'@param IEC_class a string indicating the IEC wind class (see IEC 61400-1) of the turbine to be selected. Classes \code{'I'}, \code{'II'} and \code{'III'} are suitable for sites with an annual mean wind speed of 10, 8.5 and 7.5 m/s respectively. Classes \code{'I/II'} and \code{'II/III'} indicate intermediate turbines that fit both classes. More details of the five turbines and a plot of its power curves can be found in Lledó et al. (2019). +#'@return An array with the same dimensions as wind, containing the Wind Capacity Factor (unitless). +#' +#'@examples +#'wind <- rweibull(n = 100, shape = 2, scale = 6) +#'WCF <- WindCapacityFactor(wind, IEC_class = "III") +#' +#'@export WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/III", "III")) { IEC_class <- match.arg(IEC_class) pc_files <- c( diff --git a/R/WindPowerDensity.R b/R/WindPowerDensity.R index 21ed6cb..0606b3a 100644 --- a/R/WindPowerDensity.R +++ b/R/WindPowerDensity.R @@ -1,8 +1,18 @@ - -#======================= -# Wind power density -# default: 1.225, the standard density of air at 15 ºC and 1013.25 hPa -#======================= +#'Wind power density on multidimensional array objects +#' +#'@author Llorenç Lledó, \email{llledo@bsc.es} +#'@description Wind Power Density computes the wind power that is available for extraction per square meter of swept area. +#'@description It is computed as 0.5*ro*wspd^3. As this function is non-linear, it will give inaccurate results if used with period means. +#' +#'@param wind a multidimensional array, vector or scalar with instantaneous wind speeds expressed in m/s. +#'@param ro a scalar, or alternatively a multidimensional array with the same dimensions as wind, with the air density expressed in kg/m^3. By default it takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa. +#'@return An array with the same dimensions as wind, containing Wind Power Density expressed in W/m^2. +#' +#'@examples +#'wind <- rweibull(n = 100, shape = 2, scale = 6) +#'WPD <- WindPowerDensity(wind) +#' +#'@export WindPowerDensity <- function(wind, ro = 1.225) { return(0.5 * ro * wind^3) } -- GitLab From 9ed5377e1d0a3772fbaae72feb7e65e0cb92b957 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 10:57:07 +0200 Subject: [PATCH 08/20] Generated documentation --- DESCRIPTION | 2 +- NAMESPACE | 2 +- man/WindCapacityFactor.Rd | 34 ++++++++++++++++++++++++++++++++++ man/WindPowerDensity.Rd | 29 +++++++++++++++++++++++++++++ 4 files changed, 65 insertions(+), 2 deletions(-) create mode 100644 man/WindCapacityFactor.Rd create mode 100644 man/WindPowerDensity.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 18a9d22..b383147 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,4 +19,4 @@ URL: https://earth.bsc.es/gitlab/es/csindicators/ BugReports: https://earth.bsc.es/gitlab/es/csindicators/-/issues LazyData: true Encoding: UTF-8 -RoxygenNote: 7.0.1 +RoxygenNote: 7.1.1 diff --git a/NAMESPACE b/NAMESPACE index 08ad46b..16e2e4c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,8 +19,8 @@ export(SelectPeriodOnDates) export(Threshold) export(TotalSpellTimeExceedingThreshold) export(TotalTimeExceedingThreshold) -export(WindPowerDensity) export(WindCapacityFactor) +export(WindPowerDensity) import(multiApply) importFrom(ClimProjDiags,Subset) importFrom(s2dv,Reorder) diff --git a/man/WindCapacityFactor.Rd b/man/WindCapacityFactor.Rd new file mode 100644 index 0000000..789a901 --- /dev/null +++ b/man/WindCapacityFactor.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/WindCapacityFactor.R +\name{WindCapacityFactor} +\alias{WindCapacityFactor} +\title{Wind capacity factor} +\usage{ +WindCapacityFactor(wind, IEC_class = c("I", "I/II", "II", "II/III", "III")) +} +\arguments{ +\item{wind}{a multidimensional array, vector or scalar with instantaneous wind speeds expressed in m/s.} + +\item{IEC_class}{a string indicating the IEC wind class (see IEC 61400-1) of the turbine to be selected. Classes \code{'I'}, \code{'II'} and \code{'III'} are suitable for sites with an annual mean wind speed of 10, 8.5 and 7.5 m/s respectively. Classes \code{'I/II'} and \code{'II/III'} indicate intermediate turbines that fit both classes. More details of the five turbines and a plot of its power curves can be found in Lledó et al. (2019).} +} +\value{ +An array with the same dimensions as wind, containing the Wind Capacity Factor (unitless). +} +\description{ +Wind capacity factor computes the wind power generated by a specific wind turbine model under specific wind speed conditions, and expresses it as a fraction of the rated capacity (i.e. maximum power) of the turbine. + +It is computed by means of a tabular power curve that relates wind speed to power output. The tabular values are interpolated with a linear piecewise approximating function to obtain a smooth power curve. Five different power curves that span different IEC classes can be selected (see below). +} +\examples{ +wind <- rweibull(n = 100, shape = 2, scale = 6) +WCF <- WindCapacityFactor(wind, IEC_class = "III") + +} +\references{ +Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F. J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 + +International Standard IEC 61400-1 (third ed.) (2005) +} +\author{ +Llorenç Lledó, \email{llledo@bsc.es} +} diff --git a/man/WindPowerDensity.Rd b/man/WindPowerDensity.Rd new file mode 100644 index 0000000..b0dfd77 --- /dev/null +++ b/man/WindPowerDensity.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/WindPowerDensity.R +\name{WindPowerDensity} +\alias{WindPowerDensity} +\title{Wind power density on multidimensional array objects} +\usage{ +WindPowerDensity(wind, ro = 1.225) +} +\arguments{ +\item{wind}{a multidimensional array, vector or scalar with instantaneous wind speeds expressed in m/s.} + +\item{ro}{a scalar, or alternatively a multidimensional array with the same dimensions as wind, with the air density expressed in kg/m^3. By default it takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa.} +} +\value{ +An array with the same dimensions as wind, containing Wind Power Density expressed in W/m^2. +} +\description{ +Wind Power Density computes the wind power that is available for extraction per square meter of swept area. + +It is computed as 0.5*ro*wspd^3. As this function is non-linear, it will give inaccurate results if used with period means. +} +\examples{ +wind <- rweibull(n = 100, shape = 2, scale = 6) +WPD <- WindPowerDensity(wind) + +} +\author{ +Llorenç Lledó, \email{llledo@bsc.es} +} -- GitLab From 0103f42a26e4625bf57c931d5024e31ec898c742 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 15:13:26 +0200 Subject: [PATCH 09/20] Add vignette for energy indicators --- vignettes/EnergyIndicators.md | 66 +++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 vignettes/EnergyIndicators.md diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md new file mode 100644 index 0000000..49d852e --- /dev/null +++ b/vignettes/EnergyIndicators.md @@ -0,0 +1,66 @@ +--- +title: "Energy Indicators" +author: "Llorenç Lledó, Earth Sciences department, Barcelona Supercomputing Center (BSC)" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > +--- + + + +## Introduction + +The energy sector is affected by the atmospheric ciruclation in many ways. On the one hand energy supply from renewable sources like wind, solar or hydropower relies on availability of wind, sunshine or water. On the other hand, energy demand is affected by changes in near-surface temperature. A number of indicators derived from atmospheric variables can be useful as proxies of energy production/demand. + +Currently, this package provides two indicators for wind power generation: + +--- +1. **WindPowerDensity -** Wind power that is available for extraction from the wind flow, per square meter of swept area. +2. **WindCapacityFactor -** Wind power generation of a wind turbine, normalized by the maximum power that the turbine can deliver (rated capacity). +--- + +### 1. Wind Power Density + +`WindPowerDensity` computes the kinetic energy that is available in the wind flow that traverses a unit of area swept by a wind turbine. For an instantaneous wind speed value, it is computed as: `WPD = 0.5 * ro * wspd^3` where `ro` is the air density in Kg/m^3 and `wspd` is the instantaneous wind speed at hub height in m/s. +Although wind turbines cannot extract all of the kinetic energy in the wind, and their efficiency can vary substantially at different wind speeds and among different wind turbines, this indicator provides a simple estimation of the wind resource quality. Typically, Wind Power Density is computed over a long period and its mean value is reported. + +As an example, we simulate a time series of 100 wind speed values from a Weibull distribution with scale factor of 6 and a shape factor of 2, which represent a sample of wind speed values obtained at a single location. The Weibull distribution is often assumed to fit observed wind speed values to a probability distribution function. Then, each instantaneous wind speed value is converted to its equivalent WPD. +The `mean` and `sd` of the WPD can be employed to summarize the wind resource in that location. + +```{r} +library(CSIndicators) +wind <- rweibull(n = 100, shape = 2, scale = 6) +WPD <- WindPowerDensity(wind) +mean(WPD) +sd(WPD) +``` + +If not specified, an air density of 1.225 kg/m^3 is assumed. Otherwise, the parameter `ro` can be set to a fixed value (for instance the mean air density at the site elevation could be used), or a timeseries of density values measured at each time stamp can be used to obtain more accurate results. + +```{r} +WPD <- WindPowerDensity(wind, ro = 1.15) +``` + +### 2. Wind Capacity Factor +`WindCapacityFactor` transforms wind speed values into normalized wind power generation values. The transformation is made employing manufacturer-provided power curves, for five different turbines, as described in Lledó et al. (2019). +The generation is normalized by the rated power of the turbine (i.e. the maximum power output it can achieve). This allows for comparisons between turbines of different sizes and wind farms of different installed capacities. Beware that the Capacity Factor (CF) values provided do not take into account any losses due to wakes, electricity transport, blade degradation, curtailments or maintenance shutdowns. +The function allows to choose from five different power curves that are suited for a different range of wind speed conditions. Each of the provided turbines is a representative of a IEC wind class. Generally speaking, commercially available wind turbines can be certified as IEC class I, II, III or a combination of them (I/II and II/III), according to their efficency at different wind speeds and the loads they can withstand. The basic idea is that most turbines in a same IEC class have similar power curves, and the differences of power output can be thoroughly studied with only this set of five turbines. +As the transformation of wind speed into wind power is non-linear, it is recomended to use instantaneous or 10-minutal wind speed values as input. Employing longer period means will produce inaccurate results, as far as the wind is not steady during that period. + +Following on the previous example, we will compute now the WCF that would be obtained from our sample of 100 wind speed values, when using a turbine of class IEC II: + +```{r} +WPD <- WindPowerDensity(wind, ro = 1.15) +``` + +### References +* Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F.J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 + +* International Standard IEC 61400-1 (third ed.) (2005) + + + -- GitLab From 6a56bacff40bbaa2a4d96f999ac68bfc46f2f34e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 15:52:43 +0200 Subject: [PATCH 10/20] Added figure --- vignettes/figures/WPD_histogram.png | Bin 0 -> 8379 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 vignettes/figures/WPD_histogram.png diff --git a/vignettes/figures/WPD_histogram.png b/vignettes/figures/WPD_histogram.png new file mode 100644 index 0000000000000000000000000000000000000000..940d978f658709e5a7fa6f7bbaa459a7db5d3332 GIT binary patch literal 8379 zcmch7c|4Te`@e=lg+Yrf&5%fikYt+~TZS;Q6bjjrG-Mf!wTWznsYtSw6l2Yn>=oIU z>}1P2F_y*{WB0v>r_b{|zt8ghy}sYy>-XQe&wb8y&Us(&b6xN2o(a>}(`4JvyPtu9 zflUjkZot662xVZ{v&6y(^eiToXfrT?HuZIkG=RrFd-i}pAVx+;CMKr6d-v|!w~v{b znT3Ugm6dh>{{06I95{IJAR8MSJ3BiE2L~7o=H%q$;^N}w=00@j(BZ>}j~qF2^ypC@ z9-d>zj`8yH^6~NU^YaS`2pm6t90Gwrp-@3VK_MZb6DLjx3k#n-c~V3~L{wDt)TvWq zVq&LHpFVTu%-OSN#l^)XBqSsyCC{BZCnY5%EiEl0BXj=zd0AQ63l}cP$;n;3cu`(n zUO_=YQBhGzNl95*83u#F;cx^3p`xOqs;a7{rlzj0uA!lE>Cz=lO-&>csimc*t*w3e z@?{+z9bH{rJv}`X3Z<{Ff91**0|Ns?Lqj7YBV%LZt5>gHyLQdQ#KhFp)XdE6`t|EK zZrm_8H@C2`u(Y(idGn@~m6f%%wT+F9t*xz{o!zZlx9sihZ{NP{;NWoQ&K*Zb$Gdm$ zqS0t4CnslT=X>|=xwyEvy1Kf#xw*T$V=x#G4-ZdIPcJVoEEenS?d{{^hUkPoF-0_Uu`7baYHi%=72ZV`F1qym%287x(hz%lP>CgoK2|#Kfeeq~zq} zSFc{Ze*HQnCFRYVH>s(qZ{NOs_wHRhA9D>FMe1?d|L9>+kRX{{8#Fz`)?(;Ly;}@bK`+$jFZ$KSoDK$HvCS$Hylo zCMG8*r>3TeMB>k%Kc}asXJ%$dB+~5c?A+Yk{QUgF!osg#zsO|r;^HEOLRnf`qEe~L z%gZzxZDnO;b#--ZZEbyhePd%|b8~ZRYioOZ`~1f^G6RDER!d#Q$SYxP=z#B`?ut+3 zKy?)s!JbgUUf(?_4}wF&4<9~M5q{g@(Y^|_cxcqxDFbzMk@1p2lepQZ5Bsf)itfBM zOD)oVbU!9kzsMp;HTRyMv9@CvNQ#q7KDn&S^l9ACJr582tpb^l5gJ8oAhBOQE^$0c zSk&s81A|X$xE}}+&B<~T0lmKm$FPqH!~^~xKd2ZnuBsD?x#fuI$1oW-q$!T=p`Z3$ zaN6k-xXm-im1x5XEdp!`=o$bZ0kLW*9@7jct~08RuDld$ty( zR_c4J#&T&NrmrR^q}5F7z5B{QOW+&~id~s&U}NMyx^)YqlJcPMdhaodn2 ze$?;X%pOpKh$5sk#B8L?L6zu8jhJbjLM&a^?tTYqx9KW=L@hiafZYBZ$Q&z+s%N%a z_jUQGkZRTb`p$QsFw+?q3sFpRPi*Eq=5G zg>tiu!^FM(Ww}Vn=q;+{+Cgf6O>D(@JfqmOzHy4kPxHX;(ZFRv{M~%Y>$*(M-gM^R zq|~;nl4*=3*D&2JxD9>dqYGi07xA|UtXLI+)#sIz3{S)4j@ZoO*Qu%){BU26VfOKu zmtYJ9&Knz*NN}VXL!9*Iwon|#mhF!^x8D}?<6p$38D4*$LGg&e7zR{nE4TB9dTYT` zB0H zCgRM6I`xsQ%sIxz8G9LNtDbI9=ASjiAr(tq6L%W9cHFQQ=uyp$*)Z6=>Zd~L8jDew zSOXT9SvmN$MOo0hY4e0)2D@{r$bLqd-slw zsu}z3M0{Mv*f~Gzho#*D^^6(4qG7oqtiq3e`3hIcI|sxJ8~0dbcgEq)_U8_fG`wmZ z0Mn+%YhMr0MrdrCdu3gVV3!48rEb<9ON5!wXWxr->Z6SMOA@~mXaypv)cem_S2GQ4 zsIuDPp3DXVrc_vo0M3>>DXNAhI9$Yc@n6XiB%RH|QRMvkHlE==3zx2NfVqOPq2P8&28B)YNq14Vbcm0zU zY4a?jckhQb2-JUUmo17YI+G!JV7CLQ;WzV+koGH%fiqeWzo*6dM6vs;B?3eHGevxi z^4c{N&-L!jZgA_qP#%s<6WvEd-Q==Vt*smP>&Yu_cyUaA+8%YodBi2zpKPM(%S0#k zU9+JXo%VXZ_MOpE;@5{0k{5qi3Zo8&f>5|7M*as=Qb%Dn?0)=T|)2Ycjb0 zlJdPbJtS}Yd9N{1A$ukn{{lgy5=QxKABRE&70vRT*PT3XncBJ8yMv6+;>0b<>ct7%Poj%8H`;c9eST9s#^U{)qBxjUaSK6rvdnLuHDRN`oC zmZKvx4hhA7ynV-)p}g?vAzl<9vsyg!p$!fOqKdwcOuJ zh=0o{vx{5&pG4iAZk^^pE}wRbJ#U8=rDLplmJ%F{7z?%=*5^bHzHjC`z3V4Y&D5pj zt~Kb6Z`hj3Ukl+sTE)p3F>xsZGr|zo1dr>>7j||?Sq`XNUv8LNKDK2s1vmtphlu1G zbDtnROinD)+p-wnYj^=&N-n2wK+qbO@T1goa`8)Bs~Pj`w-9viT4(qlCz~3qTRMqy z*t}s(Hv{;ULfZ;#HDt_Ww}>gBdk;A~Yj(kNZSxUt!3S=D7J2&sI)EMN4EaYh2+8a^ z^@CH7cEmr&hY{G#G-BV@AlAA44Attm#A zeeXj;xl~zfN%wm1a)!eH5s2ys>3b4(s*T7YptH=c_?iBJgW;9rQiqwiT?_^fn&3-K z70$Ws-UZYCHv!b3JmVqDmfh(dFORgv!#uSNdpvePtdVxU4Hz|uF%5zAN+CwDVbR5UbI$^um_|hQT^ceqbEaNl$2!|2g__Yryd(jO3+2 zF<{UkfDDZSp=nilIxNB_02zQ-nzU6bEr|btCNvFSf^sQkSz}lEnt~7@*$fruiz-EM zPdEBWK1n+drYrbO)p%aeJ@$-#8h~$dVTeXP?&RLG9X3BaOW8vEH>Nra!e&GMOp@ zsBtr{m+opy&*I-a2^z$eSHZ-A8R`dc$a;bLtd6IR(E$Cp#F>OXSI3$>S}&?dM5H`K z2z#-lk~7iHqWIgu?e=}<#6}0RhaO(6|K0?@K_n^K=XC;<9@GOEehIg+(Wj-S3=`kz z%2tPRSs$veeq*kP5(3NpI}`uf@IUu*Y08h}K14iFg}RJSZ=1${KIC+Hc{Pk6xH?IfX-hANTe0 zt~vr@)xHeH9LQBlU--6&9p*R**78ZB$us+tdd=v3bcM z`Nu2ic>M6% zIG$tGV%R6z8ygCo6ejGk0@_g5gOtzG&gmZ8{I+vwZYQjGl^Ux$-U zPRklDFX^0Wp~qNN1^gp>blW<8KM7F}H}Sn8|=4fTBW&JAZP+Z>{3|##}l3 zkQ0LSU#Cm&VCXs?@tT_GBf+c50N2MfGc7Lz25@|7i)sb z-F5!jm!C!IKO76$DEo7AzRzJ&AhUDT!q5g_Pz6ReD6}|G&SY>9ps_&xw_7{#ZnQE3 zu|S}{G3Jb=|IaSx3Ct;6Jy{U~q;z7e37Z{#&==okV9?iN^<_i!^9QZI0Frw6xDUQ< z0e&(U@JHmxDwOU`f|gn0$ZIB{x&>b#@a$ewBy-1`6w_i_IRtR#@70RiuCGSHqe@9{ zqCC9e0HiC4fXT-BHYh}ND;=)%r_ot0?)~R>^6ER z)8tGpgc6~#3helSFa1~=6B>LhHE@eUL0+%+VC0wvU&@1{=ulyx7WD>D#!1@SD>I}M@V4eh-hF!JV2pGvK$gX8q;B~dD5R}-;&e?Wvjak~9*7g&BQ zw%A!GPPu@p0{G72=oxxWIE<@dZz{JE!h78r_bl2~5K~RTE!ZN@b;E)!Wr-cZCBJZ) z5&~Jw^vvPg$;35scae(u&Wv zdt-8UGJhp5I@@EojD(}V+L-v!5C$jmD=Y7&u|Cs+bL6_mY5D&)um2^*|EK1!BgodH zr`=BAt$|QY&snN=PP>E#vXn)7a|9wBhY+4%ujyRB@{toM1r7yrf|oaZzZ^(xALt2Q z5sZ>Fnz{u23Md62kJ(k|_wOW~ytET6sXQ#|59*Qw?-JYwR|httwK|uVZvQb2e}&Dv znl)pyBwohz*U77eJ9ht;>~@mM{?u|Qaw4+$cgj(N@=-rUnB|fZfv`g#F~L4jIRAF? z1a;A=0MKwqzz~yUrU@JA>U8W;TeQedI@%~L19a>IaonH_>5QK9;jKdJM+o6QE+q9E zSugF#R^b?wSk)$of7QoTuTfm)!G31LK_f9OAXyF2s`!uOr3E)2cK`xAY0Hlt$Mg^} zR@=EV9*i9HTBFJDsP)9GfaGUDHfIrz1bD)RtNEJ>lc1680A*x*U#(oO64u zmb)#nA+K5U!5n>87|iOU0j)Tqu&_B3{3Y3qsS8EwREb$iDn0Rs^gAY3`tC!XhX&Q| z5ACnJQ?c6bTtt1SEG4mW{v{a1+Vy8i5Og0T@9?zqNe-(Z9+5rUUgiCzFJnIU7QK?# zv)7tF`fBGGwv#B~(D+^NpPkxvi95sI;{e@gvyu+G*L4&_kTQx1~ zq)zESN|WC;CB15C9rFjY7QG-*jIc!eDE_mq|8LOuVz;%;?(7*$bm0ezyz7thm68EH`5ba`<&V%Y=OQmfA~%9ZK0uS8D1N3_}S zuh*di9z}Ui5aI8T6qm`xZ%J0R6YQz}X36XsskAPEIK|}$h&!bvTcA0`K@82Ican#= z)wHa;?Z#cQ`f-K5cYiIP0ZK~Sm08-!T^DUL7e3P@?K!&_ch(J_+X$Cm7~1yhoxB%N zpSd*=;FKZ<8BaLF(|a&++&z;tJXm#Ec^f~;o|*!~baoZkHIZ6+e=GvZj(4owwOBy5 z2kYE%B`49&Fo6~vac^#xRNC7qD`=VPIybwbOsx1!Eb}-%TI%_&yVd6n$|GvL)U(j> zu%j!aw>6(HfmZ9=PED?9&k12K{aoaxaZlz*;_4}QEun%mjI)zazCkP9_=H8zU;I>5 zP_X{`St3Wp_N89qifI_J;usJV^<0MdNU-F1CCm9i+-^7bli;`d1bp0^;ng+M+!$y(1 z0=BrT)rQiJPh30EtV!yBcn(s+_>U@Wd`U(oUJCMh@qKuAq1W(I$>o<5ZF`f>L97s#yu_g`sxj_><{c-&iW{j;nv*#kUg|RT4VBG^R0@sVA*$-i$o( zv$DaBzmA(G)A}?+&4LPV&2dZWAI!T?Jgg3Ns}oT;>p)5jqtGJ1k9mepZ0vKjG8i<{R+nVOJ zorNFSmyd}7!ttRxwYu9C?>ld22e^mlcxePI{)9PVR(UHnW6dqkEE6WbCLuDCI^lx> zX|TDxYx_vDX*W91^&K@6?(kV=;4`gr(Yn4RZ7wIGbWO1VQ3%^0+d|Do49sT@n;t3S zbSgChN@m*>714V8Qf|$c+N4`tE%8|I2R_AyNOWlC%2!^3tcBP@6^X^`m))R}!Bo(* z&+kR?_NZRX8;Wkm0f`4~tLD7XaKAaUSuYsra$~Kc>e`Ih$beD4`L)^Jyon$E0DM%H z@&*mBS=Q8?yVVT`Bua^EMo`U@A}PtoDerq0xREYO?mgAOFFhd}Q^G?VC!gyv=O?s5 zHKg4;QUfp$iq2_#B>|CBW= Date: Fri, 23 Apr 2021 15:53:08 +0200 Subject: [PATCH 11/20] Add figures, add more content. --- vignettes/EnergyIndicators.md | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index 49d852e..3d80e1d 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -28,16 +28,22 @@ Currently, this package provides two indicators for wind power generation: `WindPowerDensity` computes the kinetic energy that is available in the wind flow that traverses a unit of area swept by a wind turbine. For an instantaneous wind speed value, it is computed as: `WPD = 0.5 * ro * wspd^3` where `ro` is the air density in Kg/m^3 and `wspd` is the instantaneous wind speed at hub height in m/s. Although wind turbines cannot extract all of the kinetic energy in the wind, and their efficiency can vary substantially at different wind speeds and among different wind turbines, this indicator provides a simple estimation of the wind resource quality. Typically, Wind Power Density is computed over a long period and its mean value is reported. -As an example, we simulate a time series of 100 wind speed values from a Weibull distribution with scale factor of 6 and a shape factor of 2, which represent a sample of wind speed values obtained at a single location. The Weibull distribution is often assumed to fit observed wind speed values to a probability distribution function. Then, each instantaneous wind speed value is converted to its equivalent WPD. -The `mean` and `sd` of the WPD can be employed to summarize the wind resource in that location. +As an example, we simulate a time series of 1000 wind speed values from a Weibull distribution with scale factor of 6 and a shape factor of 2, which represent a sample of wind speed values obtained at a single location. The Weibull distribution is often assumed to fit observed wind speed values to a probability distribution function. Then, each instantaneous wind speed value is converted to its equivalent WPD. +The `mean` and `sd` of the WPD can be employed to summarize the wind resource in that location. Otherwise, we can plot the histograms to see the full distribution of values: ```{r} library(CSIndicators) -wind <- rweibull(n = 100, shape = 2, scale = 6) +wind <- rweibull(n = 1000, shape = 2, scale = 6) WPD <- WindPowerDensity(wind) mean(WPD) sd(WPD) +par(mfrow=c(1,2)) +hist(wind, breaks=seq(0,20)) +hist(WPD, breaks=seq(0,4000,200)) ``` + + +As you can see the histogram of the WPD is highly skewed, even if the wind speed was only a little skewed! If not specified, an air density of 1.225 kg/m^3 is assumed. Otherwise, the parameter `ro` can be set to a fixed value (for instance the mean air density at the site elevation could be used), or a timeseries of density values measured at each time stamp can be used to obtain more accurate results. @@ -48,13 +54,17 @@ WPD <- WindPowerDensity(wind, ro = 1.15) ### 2. Wind Capacity Factor `WindCapacityFactor` transforms wind speed values into normalized wind power generation values. The transformation is made employing manufacturer-provided power curves, for five different turbines, as described in Lledó et al. (2019). The generation is normalized by the rated power of the turbine (i.e. the maximum power output it can achieve). This allows for comparisons between turbines of different sizes and wind farms of different installed capacities. Beware that the Capacity Factor (CF) values provided do not take into account any losses due to wakes, electricity transport, blade degradation, curtailments or maintenance shutdowns. + The function allows to choose from five different power curves that are suited for a different range of wind speed conditions. Each of the provided turbines is a representative of a IEC wind class. Generally speaking, commercially available wind turbines can be certified as IEC class I, II, III or a combination of them (I/II and II/III), according to their efficency at different wind speeds and the loads they can withstand. The basic idea is that most turbines in a same IEC class have similar power curves, and the differences of power output can be thoroughly studied with only this set of five turbines. -As the transformation of wind speed into wind power is non-linear, it is recomended to use instantaneous or 10-minutal wind speed values as input. Employing longer period means will produce inaccurate results, as far as the wind is not steady during that period. +Notice that power curves are intended to be used with 10-minutal steady wind speed values at hub height, which in modern wind turbines varies between 80 and 120m typically. As the transformation of wind speed into wind power is non-linear, it is recomended to use instantaneous or 10-minutal wind speed values as input. Employing longer period means will produce inaccurate results, as far as the wind is not steady during that period. -Following on the previous example, we will compute now the WCF that would be obtained from our sample of 100 wind speed values, when using a turbine of class IEC II: +Following on the previous example, we will compute now the CF that would be obtained from our sample of 1000 wind speed values, when using a turbine of class IEC II: ```{r} -WPD <- WindPowerDensity(wind, ro = 1.15) +WCF <- WindCapacityFactor(wind, IEC_class = "II") +par(mfrow=c(1,2)) +hist(wind, breaks=seq(0,20)) +hist(WPD, breaks=seq(0,4000,200)) ``` ### References -- GitLab From c8b91091b23cca31698df1750f2871614355a0c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 15:53:59 +0200 Subject: [PATCH 12/20] Adjust figure size --- vignettes/EnergyIndicators.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index 3d80e1d..f8496c3 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -41,7 +41,7 @@ par(mfrow=c(1,2)) hist(wind, breaks=seq(0,20)) hist(WPD, breaks=seq(0,4000,200)) ``` - + As you can see the histogram of the WPD is highly skewed, even if the wind speed was only a little skewed! -- GitLab From ce1e1af52feba3c7359cb7c18a5593ee187e0d2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 16:40:58 +0200 Subject: [PATCH 13/20] Added second figure --- vignettes/figures/WCF_histogram.png | Bin 0 -> 8993 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 vignettes/figures/WCF_histogram.png diff --git a/vignettes/figures/WCF_histogram.png b/vignettes/figures/WCF_histogram.png new file mode 100644 index 0000000000000000000000000000000000000000..37fcf4668d0fd7cff157a273fb3519546bd2502b GIT binary patch literal 8993 zcmc&)c|6qb_E*Z1kYo=bq(n&8;oD3|wr^Bq4I$ZghB2jy>`9Sq-%p4ZWQlAumXLJ@ zA?sMeN0td=8_WDY)c4-L-+O=e{_gAdyZ86k^Loyl=RD^*&pGdN&hwe@U-h+^5AYtK zp`l^czIefqhK3fTq1i)bpoJ(&E31MKf#z4;D;m&r&z?Pd_wJ>orQNr0-~Rpk>FDU_ z>FF657#JBD4;(na#Kd&);6Y|)W)>C}R#w(ShYqo^u^m2qn4O*d$dMx)92}gSoLpR7 z+}zweJUqO-ynK9oM~@!m=jRs?5D*j;Ja+7ukdTnDu&{`T$noRHMMXu$#KgqK#ZR0# zAt51g^5n@=r%p*qN}fJ_T1rYvT3Y(dnKLpnGP1I=a&mHK&z_Z+mp^yzoPvUaqN1Xb zl9IBrvWkj|s;Vjs22)d0Q&(62<(FU1pFaa}avjE#*=OiZp{ziw)3 zYG!6;ZfTOT+1c6K+dDWo+`4tk(b4ht?b}XH zPR`EGE-o&vuC8uwZtm{x9v&W^o}OM_Uf$l`K0ZFazP^5be*XUc0RaI>Br-5CFeoSp zg+c`f2cyyGJ9qBfy?Zw#B;?+`dl(Gn{{8!*p`l@6VGkZW2oDdBh=_=cjEstkijIzs ziHV7gjeYp=VO(6?qeqY8PoE|wB|UrgEIBzjB_$;_H8m|Q zEj>N`x8HsP0O0xa=NTCpnVFecSy?Y$yvWYZ&dJGn`SRtfSFc{be*Nano7~*oyu7^p z{QS3X-xd@U6c!d16%`d17r%S=uB4=-w6yg7`}bvKW##4N6%`egm6cUhRah(*hr?A@ zSJ%|k)YjJ4)zy9Y@S(oGzM-L^v9a;v$B#`-P0h{CEiEmrt*vcsZSC#tcs!m!Aarze zbar-jb#;CE^r^eMyQin8x3~B6=g)n8ef|CY0|NttgM&juL&L+vL?UryWMp)7bZl&F ze0==NmoF0&6O)sZQ&UqQ2u@E=fBpJ(W@ct~c6M%VZhn5AL?V6r_HAKdVR3PBX=#Z} zCND29udJ-BuCA`Ft*x)GQz(>;jg8IC&8@Ai?d@$Ul{)s;0!2e3IH-N${1yL{xuJkN zM<$Z)EQH$^g|iu7kG{+<(s1GI`SUf$B96w02z$9JM=8eSWGwh{Tz+)I(){`J^Kl<& zGcMRvYFN|e2*qjebFdc)Jjgj8&f@f-;X#Pez);pF1?KMMYrSSx#zuIz)!6OMt1zo7 zvPc&}*I1FJ#GfN;nC8wJoASfG7?w;LfL7uc@RtLJ{^_*oz*byvmP(*b-PvMiws_r( zpx?DOKBgzP*P7er3w?UDYS}l?T@_pPeb=RpD|1T9WZ7L3P_I$~PYL+2@Qp|+y|3F| zRq;9{+g|q424Oy+_x7pS^fuZ&c_r=*!NFx_)8%uFmReD?XU zz6@|FwS@DhGL`ZWdq$u6mGGY>6-!sCmW}}l38jYo5f$W3k@4Lel z&>~(ddK?PolnH7wc8;W#sYE_e6oj?Wm}sJvBObCaE2!=J7D-JtTL25~W;%J1lJUBJ z;Pt+^_Lrb)u*)PaReaU6$8~vNDHv3J$TMxM=G#BdyU=ajg&ozU#K>)jCCLqL43P+= z!HKW0<{en@q1e==OZA0pw5b7(EU;$D0c=1oON(zx)oZW$yfUp8RNcC#DI(VH(KHJx zw@6aKa0)kr>z9`6j7KU7I>V--@{-oyR5>=bD!>QQJgeV25=c-HU-&X6A@ZyH_({)>Cy^oDl`caI$NKDPf~ zo~9F~dRg-FQXTwG(Tt6;BZuf39i56RIA}Z;7*QXe_sB=`&ARUY#XJb^mOAJ?i@8=| zw&xI6Izie!k2JjD`(sA0bT!5(zG8j4c0Uv6;8JKk%*AmnUQgp-k)j&l7jWZqxx6b# zO;L61%&h4;J@2Sfm+pMh=4t%`M!x)x{E8KPnq~S0|IEe5)cFi~gidKWP z!Mx*S-+#|}|9`FhcX$7Q1o_clP-%;?uj<_VnSR;1e~31qJ~b4acXBhdl;-bWc`R|W zzoz@$Ql&hg^-x_Go+w0K>vnNZDbaM(`>U@ZIN#-69~{YLgoSKcA1*%#H0Jww8jXLr z_A^tZrtt?S^4KBdW^3XTTd1_$n8htf!C705P7b!7Iprtp9sVeExUj{&uJ&7nH0uxc z$7jjefvVF3=4y*Dz{0Uw5B5Rp?G5cH?FZV_90`(a&2Y51bO-0g^Fx z%0rD@Cr$-_C%Ey2T^AE4?_#OhUwCSrexS@goF?(xTF8E;b}VcSV6b?w15FpdIgweV zcvXC7<9g}Uj{r=&z%(lV1CjjeqP{~@5Io)Vdg5pVCHkxqlSD_ zz`bZ%+;ZvcKStB>M5-hAuDeDP99R`nj6tQ54~?GIm^kywS;9W4~nvJ8}<3Kpo0bV-+OqMnav!um)(I zR2Kx3Ypqxn@n(cl$$H@;kR<0|+2t^9fe}bhYUPrvODSUobf9qckIsA03)VxC4?7D| zG6IJte+qQcXMD`2gd*3?>}3)0E*OE}$>bXZ2x43m6}ec)%{Kb`f2QLEU`shv-XMd!bqp4g z{!n)1RD}?u>AqS}Lvl<}u@}_W<~`pRli6)ATY^JF^WH!$_P94^Sgf^wHgqsG<6$)4wJ?KDfGKQ zR_?5EL5fAgEtNyaYdc5#v40h^;Zc^x)=L$@?c`sCdEk6-w`)JW+|T+aW8||ap*ySo zOh%@9VZDNC{$K|W;+dtvKKv&R)hI%J2vlF#)IgVr+9pul{U!!MB3}xqm~Yw&Vnhq@Dhw(v`gyd?X*^ z4}nAPzFGRYKvOo7Pj4+WU9j}ifTm2)v>|Ya1Dc-hBPg^oK(@%*2bn>lX2-WyRVX<# z3B|=n{%oSM^d9p+9-|N-eytSn<_2^aL$yJ}kR3amBO}E}EJ86-;@o22>9Ux-C;|st z`7>v5=mo!`OMmf>jUL~Cd{Y$SX}O6=Pd0)CCmuc+|6P%!h!I2~!w-1;&twP)DkgVp z`l+A5Am8=gc(FY4r_{_>nB&tnM)z}RkQ<-OD_-(5)3|)Ouw@2U*mn%8)RfkTBS?Pr zx_CrHui1ggmS!2~aRs$P^ydP?>{z^9?VmnTVJp7?M-XiXyVKHkm~!JG;1WYW-Hefn z|3nZR3JGq@BUx;f$hH0iOjh@n^yeHfc_KfG=+^<|bHiSAKdTbB6~daTrBD8ADe)gF zo4;dPI$;UjBledJi`>Td|DM@`;CQO)!9S$hMtyg5kNvN0E~X~0W?>R}3XAS5S+KyHBR-cHtiU;y4rIobGfp(ki2Z#NG^6bk|F5?5xR@pldb ze9~LWWEKSBkRL+b!YH%cztlEdRKj2pS=p;T5nNWnyE-)skx%E+{o|zgA2s#ghRG7+ z)>l{{5=%8_R4}==+o?dL2==J@*z^;lbS($|4c*c|Ol5FzA=ygheYtO1so3w`odB+! zIi{-z{3+_(zIV_+SMBfgU2ms+{`wPBY%dt$A8iS>kWl*focvs1i~m~oqauIQdZy7= z_(y5r8K9%n^879zNtpb$)jF@~V9~@ec;Yeg+5)cyTYg6*;`auUyr%b~SN>6SkM}~Y zvgu)p&<2m-Kxf}#XuKr)sUyo)mVQJ0&3uOg0y-e)o7=d@nBlE+UNi}-$+rUPbiKgS zjeJ}|{kHpKPnQz%2XZ^TK1Nz=#WU=9P3cEMzOyW+da&x3n~&}iZP z!=yQ5@|wiSRbc_I+D}ft%BhP5^MtV(hAd3qh zr4iiZV-Kjjcq9asz4OB44}TatS+<6>L!r>yRoJ2*P`|vTgBQS9?%VNosbvD%S=FBl zF1>fBLgWI(S^k?*A0*(0oZlB=R1AdSer(7w3`@N2z~9B_gRdaznS{0`sbk)7V~ro28C;Z)RZ+OxGb{OrEsdrR z?U+$Zn*1!h;W?G4(cf6@WMq}Mh%q5TP>Ki{t@RQrtoTpDr90`Ga zcR62hk``b7ejnGYwgX61B%G1yq_KrO?-CCs6V`^dF78~(%9it9YEicobi5-IDM3yR z$rb_)=>9ODyTN#SPK%zeQ@*pZ(wBflMV?JJr2_PwZ7KafScN(NV z22>oLwgpbzcrY3|wJeoNU`!=M4~Q~sF^?y2h0&=^E3FDt>4BW=$e8{_8=-+#ojR#i|N1nKA;j=pnBN(6Nwr$T*h9bQ7Ey0_;{D+CfGfpJ0mQYl@%tmy znAC(=CO9eb4K^DNRRgDC(JK1o)g}imT+O8bc4RT3*YZ|y!-&}xTglLUvMEZ{eF>C`TI${?L=vdFVTG?WCX5S37Wkdw-rLfBF-?`X9 zi3S^74lo7P>gvYkpby{C?A(s3D$n`|j>NX@u{7}*`Qg5zRiB2jnT9!utI$=XKXao) zey%`K_cZX=8OU%}4BeJ2=PkwaL4)qf`l13^hEpp_hLfLkj&p;jR9YQT&*FF6&PdRgi>fban2V(6i4;g& z4_a3Vz9~tUYA5N7?J}8l8&q@2EC5y8OuVBNp;uO%aF{=koY>R$eyixVZO4{py3WQ` z16vPr-(*(XN?EZHUsV0@WL6Su&=a0IAea;tNElpTkPCGivXN|cXniWwU~}a%zGcHNNtIcS)~WPmNXe z_MHOm&m&7F85?|&HHfxG<)*2}ardG$x)&O1gWgPTD5&R$L|H2v$;NuvEHzXF+l+U9 z9k4LX2#D~gC}=jh(#sdA)-V^}Y82Yd6IrLV5wsq9x8eDaoW1!6Yqf-paz1nk1=WBk z9XMa1!*_+`bp;xfKbw$Q(T;OAii0M_k$_tq48fB;tI-POlE1VTo%*W&RpFC6F?t z`O)a8Lwk{s!FrD4mwZH8A@Yjt1Keiors)N=^2FS+NymhxGj*YA0Ue?J53`HqrtUI` z<6ym!?Yg7iu1_}XEjWETE=#r6pZAK}*d$*hFUSw~Mop(TfG|EGR zMJE>|`m7>BgX~KF52|gR0`{K?$)>3_LI#cwwclQyaEs50e1Et61mD_Q7ayNbGpp}o zk9=JiyjL{95YSBNFm6guo|h2Euk}Sd%&zP}RY#8*72-HM()9yW9t1p_P;Gq+DCBv3 z!56lfRJBiIfwis~Vb|ZJI_c(LTW~B>VDlpPk(GGXc4v_<62@vO7Tbf~29EFTCpv@- z%A9)^>*}-L|CW96w8@JCwLIC>}B}_zMROyyIS~l$;IWsL9GMJ z{_{0`E*RB`N0(ugNa7b#n&*O1yY(3LSb27o`83&5nrNEJzBT>0^xo6#>+N+1x57M_ zqnF@mo$b^sy{!@MNQ)25pEUZa!)e2qBavHavm&ryqnt)BgC21)rsq5nUxQ}{j8GG*XIY4mk3UT#A=BWg|+`M!O{3;^f zkd&RGhb7&^PAP}AQNzbmC^dtQNmS9h_}0a@WtR^lk7Q$2MrtK=WC6VUZA~ewQmvP( zM(qhP0c$J$v9Y}I8^6|;L1}9lxi>v&!rCk2IQfJ<{{0(Zp-UIkNuh}FJN&Yb0d>*|8d7=e0@Y$XzU&Z@m0@hV#pBhX~8sSVOxy7xc*DtMPl auT){!>OOnj9{R5fjkbpVg*^C;kpBY8l6ai} literal 0 HcmV?d00001 -- GitLab From d5ed8b2a3f5076d135d31cc1e9bf207d3f00575e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 16:41:45 +0200 Subject: [PATCH 14/20] Added second figure --- vignettes/EnergyIndicators.md | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index f8496c3..5eb8597 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -37,9 +37,9 @@ wind <- rweibull(n = 1000, shape = 2, scale = 6) WPD <- WindPowerDensity(wind) mean(WPD) sd(WPD) -par(mfrow=c(1,2)) -hist(wind, breaks=seq(0,20)) -hist(WPD, breaks=seq(0,4000,200)) +par(mfrow=c(1, 2)) +hist(wind, breaks = seq(0, 20)) +hist(WPD, breaks = seq(0, 4000, 200)) ``` @@ -56,17 +56,24 @@ WPD <- WindPowerDensity(wind, ro = 1.15) The generation is normalized by the rated power of the turbine (i.e. the maximum power output it can achieve). This allows for comparisons between turbines of different sizes and wind farms of different installed capacities. Beware that the Capacity Factor (CF) values provided do not take into account any losses due to wakes, electricity transport, blade degradation, curtailments or maintenance shutdowns. The function allows to choose from five different power curves that are suited for a different range of wind speed conditions. Each of the provided turbines is a representative of a IEC wind class. Generally speaking, commercially available wind turbines can be certified as IEC class I, II, III or a combination of them (I/II and II/III), according to their efficency at different wind speeds and the loads they can withstand. The basic idea is that most turbines in a same IEC class have similar power curves, and the differences of power output can be thoroughly studied with only this set of five turbines. + Notice that power curves are intended to be used with 10-minutal steady wind speed values at hub height, which in modern wind turbines varies between 80 and 120m typically. As the transformation of wind speed into wind power is non-linear, it is recomended to use instantaneous or 10-minutal wind speed values as input. Employing longer period means will produce inaccurate results, as far as the wind is not steady during that period. -Following on the previous example, we will compute now the CF that would be obtained from our sample of 1000 wind speed values, when using a turbine of class IEC II: +Following on the previous example, we will compute now the CF that would be obtained from our sample of 1000 wind speed values when using a turbine of class IEC I, and compare it to the CF values for a class III: ```{r} -WCF <- WindCapacityFactor(wind, IEC_class = "II") -par(mfrow=c(1,2)) -hist(wind, breaks=seq(0,20)) -hist(WPD, breaks=seq(0,4000,200)) +WCFI <- WindCapacityFactor(wind, IEC_class = "I") +WCFIII <- WindCapacityFactor(wind, IEC_class = "III") +par(mfrow=c(1, 3)) +hist(wind, breaks = seq(0, 20)) +hist(WCFI, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) +hist(WCFIII, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) ``` + + +Blabla. + ### References * Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F.J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 -- GitLab From a821545ec46d115f361635e2138e83dff7821131 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 16:42:23 +0200 Subject: [PATCH 15/20] Adjust figure size --- vignettes/EnergyIndicators.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index 5eb8597..799babb 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -70,7 +70,7 @@ hist(WCFI, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) hist(WCFIII, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) ``` - + Blabla. -- GitLab From 4608db5980e21c74094f44c4a7a845ef4cbe60f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 16:43:17 +0200 Subject: [PATCH 16/20] Forgot to change figure name --- vignettes/EnergyIndicators.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index 799babb..60dbe99 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -70,7 +70,7 @@ hist(WCFI, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) hist(WCFIII, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) ``` - + Blabla. -- GitLab From c25bacc288f5b7635fef29732402ed598be5fd51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 17:00:03 +0200 Subject: [PATCH 17/20] More content --- vignettes/EnergyIndicators.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index 60dbe99..cd8af47 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -72,7 +72,7 @@ hist(WCFIII, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) -Blabla. +From the CF histograms we can see that, for this particular wind speed distribution, the IEC I turbine (designed for high winds) producess less energy than the IEC III turbine, which in is more suitable for this range of wind speed values. ### References * Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F.J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 -- GitLab From 1654742907c2ed5a4fbd311c217b76fe13f60ad3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 17:04:04 +0200 Subject: [PATCH 18/20] Correction --- vignettes/EnergyIndicators.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index cd8af47..88dfa49 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -52,7 +52,7 @@ WPD <- WindPowerDensity(wind, ro = 1.15) ``` ### 2. Wind Capacity Factor -`WindCapacityFactor` transforms wind speed values into normalized wind power generation values. The transformation is made employing manufacturer-provided power curves, for five different turbines, as described in Lledó et al. (2019). +`WindCapacityFactor` transforms wind speed values into normalized wind power values. The transformation is made employing manufacturer-provided power curves, for five different turbines, as described in Lledó et al. (2019). The generation is normalized by the rated power of the turbine (i.e. the maximum power output it can achieve). This allows for comparisons between turbines of different sizes and wind farms of different installed capacities. Beware that the Capacity Factor (CF) values provided do not take into account any losses due to wakes, electricity transport, blade degradation, curtailments or maintenance shutdowns. The function allows to choose from five different power curves that are suited for a different range of wind speed conditions. Each of the provided turbines is a representative of a IEC wind class. Generally speaking, commercially available wind turbines can be certified as IEC class I, II, III or a combination of them (I/II and II/III), according to their efficency at different wind speeds and the loads they can withstand. The basic idea is that most turbines in a same IEC class have similar power curves, and the differences of power output can be thoroughly studied with only this set of five turbines. -- GitLab From 0beb684088a316398d4da0db9d666a1fc16e5141 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 23 Apr 2021 17:06:06 +0200 Subject: [PATCH 19/20] Correction --- vignettes/EnergyIndicators.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.md index 88dfa49..c02ec8f 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.md @@ -72,7 +72,7 @@ hist(WCFIII, breaks = seq(0, 1, 0.05), ylim = c(0, 500)) -From the CF histograms we can see that, for this particular wind speed distribution, the IEC I turbine (designed for high winds) producess less energy than the IEC III turbine, which in is more suitable for this range of wind speed values. +From the CF histograms we can see that, for this particular wind speed distribution, the IEC I turbine (designed for high winds) producess less energy than the IEC III turbine, which is more suitable for this range of wind speed values. ### References * Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F.J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 -- GitLab From 280078e0e2861e1639d2cac323c9d40ff8581cc7 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 29 Apr 2021 09:48:27 +0200 Subject: [PATCH 20/20] CST version figures compression --- DESCRIPTION | 11 +++-- NAMESPACE | 4 ++ R/WindCapacityFactor.R | 37 +++++++++++++++++ R/WindPowerDensity.R | 34 +++++++++++++++ man/CST_WindCapacityFactor.Rd | 39 ++++++++++++++++++ man/CST_WindPowerDensity.Rd | 34 +++++++++++++++ ...ergyIndicators.md => EnergyIndicators.Rmd} | 10 ++--- vignettes/figures/WCF_histogram.png | Bin 8993 -> 5814 bytes vignettes/figures/WPD_histogram.png | Bin 8379 -> 6445 bytes 9 files changed, 161 insertions(+), 8 deletions(-) create mode 100644 man/CST_WindCapacityFactor.Rd create mode 100644 man/CST_WindPowerDensity.Rd rename vignettes/{EnergyIndicators.md => EnergyIndicators.Rmd} (97%) diff --git a/DESCRIPTION b/DESCRIPTION index b383147..dfea00e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,8 @@ Package: CSIndicators Title: Sectorial Indicators for Climate Services from Sub-Seasonal Forecast to Decadal Predictions Version: 0.0.1 Authors@R: c( - person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071"))) + person("Núria", "Pérez-Zanón", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")), + person("Llorenç", "Lledó", , "llorenc.lledo@bsc.es", role = "aut")) Description: The package contains the definition-based computation of the sectoral indicators for the Climate Service. Depends: R (>= 3.6.0) @@ -13,10 +14,14 @@ Imports: ClimProjDiags Suggests: testthat, - CSTools + CSTools, + knitr, + markdown, + rmarkdown +VignetteBuilder: knitr License: Apache License 2.0 URL: https://earth.bsc.es/gitlab/es/csindicators/ BugReports: https://earth.bsc.es/gitlab/es/csindicators/-/issues LazyData: true Encoding: UTF-8 -RoxygenNote: 7.1.1 +RoxygenNote: 7.0.1 diff --git a/NAMESPACE b/NAMESPACE index 16e2e4c..a4b53eb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -11,6 +11,8 @@ export(CST_SelectPeriodOnData) export(CST_Threshold) export(CST_TotalSpellTimeExceedingThreshold) export(CST_TotalTimeExceedingThreshold) +export(CST_WindCapacityFactor) +export(CST_WindPowerDensity) export(PeriodAccumulation) export(PeriodMean) export(QThreshold) @@ -24,5 +26,7 @@ export(WindPowerDensity) import(multiApply) importFrom(ClimProjDiags,Subset) importFrom(s2dv,Reorder) +importFrom(stats,approxfun) importFrom(stats,ecdf) importFrom(stats,quantile) +importFrom(utils,read.delim) diff --git a/R/WindCapacityFactor.R b/R/WindCapacityFactor.R index dba82f9..2b60ae4 100644 --- a/R/WindCapacityFactor.R +++ b/R/WindCapacityFactor.R @@ -1,3 +1,38 @@ +#'Wind capacity factor on s2dv_cube objects +#' +#'@author Llorenç Lledó, \email{llledo@bsc.es} +#'@description Wind capacity factor computes the wind power generated by a specific wind turbine model under specific wind speed conditions, and expresses it as a fraction of the rated capacity (i.e. maximum power) of the turbine. +#'@description It is computed by means of a tabular power curve that relates wind speed to power output. The tabular values are interpolated with a linear piecewise approximating function to obtain a smooth power curve. Five different power curves that span different IEC classes can be selected (see below). +#'@references Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F. J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 +#'@references International Standard IEC 61400-1 (third ed.) (2005) +#' +#'@param wind a s2dv_cube object with instantaneous wind speeds expressed in m/s. +#'@param IEC_class a string indicating the IEC wind class (see IEC 61400-1) of the turbine to be selected. Classes \code{'I'}, \code{'II'} and \code{'III'} are suitable for sites with an annual mean wind speed of 10, 8.5 and 7.5 m/s respectively. Classes \code{'I/II'} and \code{'II/III'} indicate intermediate turbines that fit both classes. More details of the five turbines and a plot of its power curves can be found in Lledó et al. (2019). +#'@return A s2dv_cube object containing the Wind Capacity Factor (unitless). +#' +#'@examples +#'wind <- array(rweibull(n = 100, shape = 2, scale = 6), c(member = 10, lat = 2, lon = 5)) +#'wind <- CSTools::s2dv_cube(data = wind, lat = c(40, 41), lon = 1:5, +#' Variable = list(varName = 'sfcWind', level = 'Surface'), +#' Datasets = 'synthetic', when = Sys.time(), +#' Dates = list(start = '1990-01-01 00:00:00', end = '1990-01-01 00:00:00'), +#' source_file = NA) +#'WCF <- CST_WindCapacityFactor(wind, IEC_class = "III") +#' +#'@export +CST_WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/III", "III")) { + if (!inherits(wind, 's2dv_cube')) { + stop("Parameter 'wind' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + wind$data <- WindCapacityFactor(wind$data, IEC_class = IEC_class) + if ('Variable' %in% names(wind)) { + if ('varName' %in% names(wind$Variable)) { + wind$Variable$varName <- 'WindCapacityFactor' + } + } + return(wind) +} #'Wind capacity factor #' #'@author Llorenç Lledó, \email{llledo@bsc.es} @@ -10,6 +45,8 @@ #'@param IEC_class a string indicating the IEC wind class (see IEC 61400-1) of the turbine to be selected. Classes \code{'I'}, \code{'II'} and \code{'III'} are suitable for sites with an annual mean wind speed of 10, 8.5 and 7.5 m/s respectively. Classes \code{'I/II'} and \code{'II/III'} indicate intermediate turbines that fit both classes. More details of the five turbines and a plot of its power curves can be found in Lledó et al. (2019). #'@return An array with the same dimensions as wind, containing the Wind Capacity Factor (unitless). #' +#'@importFrom stats approxfun +#'@importFrom utils read.delim #'@examples #'wind <- rweibull(n = 100, shape = 2, scale = 6) #'WCF <- WindCapacityFactor(wind, IEC_class = "III") diff --git a/R/WindPowerDensity.R b/R/WindPowerDensity.R index 0606b3a..a231e45 100644 --- a/R/WindPowerDensity.R +++ b/R/WindPowerDensity.R @@ -1,3 +1,37 @@ +#'Wind power density on s2dv_cube objects +#' +#'@author Llorenç Lledó, \email{llledo@bsc.es} +#'@description Wind Power Density computes the wind power that is available for extraction per square meter of swept area. +#'@description It is computed as 0.5*ro*wspd^3. As this function is non-linear, it will give inaccurate results if used with period means. +#' +#'@param wind a s2dv_cube object with instantaneous wind speeds expressed in m/s obtained from CST_Load or s2dv_cube functions from CSTools pacakge +#'@param ro a scalar, or alternatively a multidimensional array with the same dimensions as wind, with the air density expressed in kg/m^3. By default it takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa. +#'@return A s2dv_cube object containing Wind Power Density expressed in W/m^2. +#' +#'@examples +#'wind <- array(rweibull(n = 100, shape = 2, scale = 6), c(member = 10, lat = 2, lon = 5)) +#'wind <- CSTools::s2dv_cube(data = wind, lat = c(40, 41), lon = 1:5, +#' Variable = list(varName = 'sfcWind', level = 'Surface'), +#' Datasets = 'synthetic', when = Sys.time(), +#' Dates = list(start = '1990-01-01 00:00:00', end = '1990-01-01 00:00:00'), +#' source_file = NA) +#'WPD <- CST_WindPowerDensity(wind) +#' +#'@export +CST_WindPowerDensity <- function(wind, ro = 1.225) { + if (!inherits(wind, 's2dv_cube')) { + stop("Parameter 'wind' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + wind$data <- WindPowerDensity(wind$data, ro = ro) + if ('Variable' %in% names(wind)) { + if ('varName' %in% names(wind$Variable)) { + wind$Variable$varName <- 'WindPowerDensity' + } + } + return(wind) +} + #'Wind power density on multidimensional array objects #' #'@author Llorenç Lledó, \email{llledo@bsc.es} diff --git a/man/CST_WindCapacityFactor.Rd b/man/CST_WindCapacityFactor.Rd new file mode 100644 index 0000000..0c75557 --- /dev/null +++ b/man/CST_WindCapacityFactor.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/WindCapacityFactor.R +\name{CST_WindCapacityFactor} +\alias{CST_WindCapacityFactor} +\title{Wind capacity factor on s2dv_cube objects} +\usage{ +CST_WindCapacityFactor(wind, IEC_class = c("I", "I/II", "II", "II/III", "III")) +} +\arguments{ +\item{wind}{a s2dv_cube object with instantaneous wind speeds expressed in m/s.} + +\item{IEC_class}{a string indicating the IEC wind class (see IEC 61400-1) of the turbine to be selected. Classes \code{'I'}, \code{'II'} and \code{'III'} are suitable for sites with an annual mean wind speed of 10, 8.5 and 7.5 m/s respectively. Classes \code{'I/II'} and \code{'II/III'} indicate intermediate turbines that fit both classes. More details of the five turbines and a plot of its power curves can be found in Lledó et al. (2019).} +} +\value{ +A s2dv_cube object containing the Wind Capacity Factor (unitless). +} +\description{ +Wind capacity factor computes the wind power generated by a specific wind turbine model under specific wind speed conditions, and expresses it as a fraction of the rated capacity (i.e. maximum power) of the turbine. + +It is computed by means of a tabular power curve that relates wind speed to power output. The tabular values are interpolated with a linear piecewise approximating function to obtain a smooth power curve. Five different power curves that span different IEC classes can be selected (see below). +} +\examples{ +wind <- array(rweibull(n = 100, shape = 2, scale = 6), c(member = 10, lat = 2, lon = 5)) +wind <- CSTools::s2dv_cube(data = wind, lat = c(40, 41), lon = 1:5, + Variable = list(varName = 'sfcWind', level = 'Surface'), + Datasets = 'synthetic', when = Sys.time(), + Dates = list(start = '1990-01-01 00:00:00', end = '1990-01-01 00:00:00'), + source_file = NA) +WCF <- CST_WindCapacityFactor(wind, IEC_class = "III") + +} +\references{ +Lledó, Ll., Torralba, V., Soret, A., Ramon, J., & Doblas-Reyes, F. J. (2019). Seasonal forecasts of wind power generation. Renewable Energy, 143, 91–100. https://doi.org/10.1016/j.renene.2019.04.135 + +International Standard IEC 61400-1 (third ed.) (2005) +} +\author{ +Llorenç Lledó, \email{llledo@bsc.es} +} diff --git a/man/CST_WindPowerDensity.Rd b/man/CST_WindPowerDensity.Rd new file mode 100644 index 0000000..a0d552c --- /dev/null +++ b/man/CST_WindPowerDensity.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/WindPowerDensity.R +\name{CST_WindPowerDensity} +\alias{CST_WindPowerDensity} +\title{Wind power density on s2dv_cube objects} +\usage{ +CST_WindPowerDensity(wind, ro = 1.225) +} +\arguments{ +\item{wind}{a s2dv_cube object with instantaneous wind speeds expressed in m/s obtained from CST_Load or s2dv_cube functions from CSTools pacakge} + +\item{ro}{a scalar, or alternatively a multidimensional array with the same dimensions as wind, with the air density expressed in kg/m^3. By default it takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa.} +} +\value{ +A s2dv_cube object containing Wind Power Density expressed in W/m^2. +} +\description{ +Wind Power Density computes the wind power that is available for extraction per square meter of swept area. + +It is computed as 0.5*ro*wspd^3. As this function is non-linear, it will give inaccurate results if used with period means. +} +\examples{ +wind <- array(rweibull(n = 100, shape = 2, scale = 6), c(member = 10, lat = 2, lon = 5)) +wind <- CSTools::s2dv_cube(data = wind, lat = c(40, 41), lon = 1:5, + Variable = list(varName = 'sfcWind', level = 'Surface'), + Datasets = 'synthetic', when = Sys.time(), + Dates = list(start = '1990-01-01 00:00:00', end = '1990-01-01 00:00:00'), + source_file = NA) +WPD <- CST_WindPowerDensity(wind) + +} +\author{ +Llorenç Lledó, \email{llledo@bsc.es} +} diff --git a/vignettes/EnergyIndicators.md b/vignettes/EnergyIndicators.Rmd similarity index 97% rename from vignettes/EnergyIndicators.md rename to vignettes/EnergyIndicators.Rmd index c02ec8f..1209c2d 100644 --- a/vignettes/EnergyIndicators.md +++ b/vignettes/EnergyIndicators.Rmd @@ -4,13 +4,13 @@ author: "Llorenç Lledó, Earth Sciences department, Barcelona Supercomputing Ce date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Energy Indicators} + %\usepackage[utf8]{inputenc} --- - +Energy Indicators +----------------------- ## Introduction diff --git a/vignettes/figures/WCF_histogram.png b/vignettes/figures/WCF_histogram.png index 37fcf4668d0fd7cff157a273fb3519546bd2502b..493f2256b3342322aaf1680e38732f9a88fc1962 100644 GIT binary patch literal 5814 zcma)A3p|u*`*sMc!lra^n%cw`k))iOK{9C=Ws^=uIj&Pvp~)~!s%4!w)*-?miP9t@ zG;GU|Gz>!0m@miIX(qfm3^Vh+??l`F_wfDw{2tHyJkNdK*L7dl^UnN49>JpJRw}O) z6BCoOIb`J^CMF4piHU>XNJ1;OXJ!%5pO7QA#|{dGLV0<4Lqo%J=gxV0dX|@$TUc1M zx3@buI2afh#K*^5Sy>qy8|UZe^ZERsprEFvrsU-0tgNh$A3xr|e?L4t+|SRixw*Nz zx|%>BP^nZFi*@qkNdy9c!C=bD%EZOR>2$ibw)Xh=_^n&F$Yk=+(9ow(pVHFO_#2Lt ziis)B*jQN{yYgyJXr+o(wA*pPQ{gknU%c z`_E_WJk9#?77AVDSjwTN9|Z4gutNv46S9*R=K_+tJcU%7hKr?`KF);8+t*d)ck#;0h1rz(=W`1~LFqNCr}K@J>E4U$CfbcEP64R*mlV&AC^tf!3-AH10@qncc)+IBgxP5&h)unxMWu zlgx)%(;n*%c7<4@Y7@GyS_@Ob?JT>JbG_=rW7uld;pqXlw$hGftQLPx;F9P$f9VIa z0X#sRG@)?}cnW3DTN1x-5ld2}A~iF- z(vM%>J<>2MOxqK?K+kev5{mkcasYPErF8ooo9ft{Pj=5Yp*g^dl;@s+@J(oCNUXmr z2Gdy;k)PiAesx*t-M=%{XPbHWXg2QferHh_WP^(0$(bJd6pmz{?lY75@MZk9V#A&G#s)_QRobip;*#ng2JtoiqeiMR^CqwwUQHyO zzmTkt($z^{yz6WZ{Me zUI~e_t$Qz8NU0e3gvs+fz_y~ITAXRH3G zDFVU#6Rtltmu7HD_=JehKbM3V{|hM*$Il&sc@(M>LjxAg2GfAN?5rx!RKBg@N@ru}dlh(GBSXduDYal0Ss`CApadWMUz&@gSZ=V4 zND)~0Wr@Db2!A85%gszM&L27s5-sr=XG-0^iu!+~{KJEA$9Y+4Hct+x!9FO$M1A7P zLCI*n%;BC16DbJA7-=U~Iks^}X;D|AjU0dw>E3%mafR;Jc%0traa;+QIg#$rF&N^|H6RV(=rSMp6cbN=s;6Xmh;Dz|Rg?Jh7O19FC}g`~auffu6z^m`R0!kU9hvh& zdC8j#L`6XSpI-=nHiyUyf4*Qmsje%%H*l;}*ma=dL~mPtza*d(3sbYuYZU}nDr>A_ zK7c$75yR!UAA2^+zJX&b9cQ3h-6vvj9ZQcz>N%v9L`S#l>KAO8O=cE8vne~750sZ; z2PXgzuK})VEMD(ge00ZU;T3zF>Ze4_hvws8sck)vCP=v@2&5;J5_lsNQK0TA(b^@H11) zuj6bXm5l_FAg_qp=Cf&MRc7TfoYxq2Z|xv!s-G?V4kt!F@1nJRuI_=g>>a}mex|)! zYYR7sdg(8RT&O7CMEo_lr=sw-lusAFSw@Tj1wCmNR3FhEuS*48;~@oO5OK>y@CD_g zA~i%G9+Vr85;-OWwHi)E<>wty5&b)owr>E6)~yB1xtL>TPlm~X$D?EeDYKqE@0IBK;zVuw z3KM3Cfu@4Lsm`_D!}?I}dLX7%mbod{z*>tCj+8?cT1hG&*EMTNA2P?^O_gbxelE*I zY$R4@IZbZLU0V;--n_bz={3Xg%FFwl0j$<2KCC20%W#!Y14uSyt8@INm&$IN?DJEb z>-@Y$<0EAOIcX6mS(w$=*uSD7TF(&sgZ$TgVd@G=+^%;Jk1*tKp*5YFP%3bXrSLLp zHM}WlyXIXZ+t*(j-!yosMjAf#Oy}Y3ZQ)NLB|)kUS#VPkXLdJLEk-yKD3YX?ja!3s z;chFKHg!E-G=eyoBOL#&{ZNJ**N~Hib6j&RbA$$%?JXt}38wRJu>Y(jtn!kuTJVmj zBVP%-JYDV-N)-$qb!ue41+Kw#75pdbx{F)IA+%69um1oW2A-b_^(rF{&I^wcfhM9V ze?MZ;aA7qlOf;?l0s|>)Dd}~NJATkrqa}&dj9jUs{di=_QH-(*?I<@(G(IoYT!x6R zAISEPkXkw!(+Zh`!uB{ZU(Du~Bh_fPT-8z0750!e0y*ve!AliQ01&&bWKL&jN)|Bp zifP`MA;!n`>~I66B_V@q1lW_p<3x^x{B(q;#%dY#KmAtBH`>?}XqoCojwk}N*a8Jo zS+dk`1b_X8)uir4Rk$dXb>$86v|g|(Gi6;KkNYoTHNpoV@u3QK(o zwoYC5hg|myQ|o+E%8DKc;;!uI22)5c{eO^(BFW!E@iMqru}KO_2SxHb&dMKjjae4? zYr&UI1H&<}hTk1nG&J_yP#^c`1>lB!n~nFRI+h_|Xi%>FK81XGGmF2l5B4z93WzK} zp;JZ{GB6ZLnj0V!yy`6X8!(xYH4wbZ?Pl{vel^p@pyEaabgxE$12@$n{{@OnhDGP+ zO7d6o>ZrRV;XTQpbjdAH&JkE^c`zghv=VggHHF&&Zn+HVth2wTLP<9#KmfVk5Jt6( zI0!0n6|?21Lt=rH_Es6${Yo`22e-f!5Ca>Xo3D^uV|tCjMq3S4oj?;_hAFQDeK{!n zs4beif>(Q|soMcT)WKr$emHV@tYycpa3slTWBcHN%bQ>^DH|Ca0(=JDE@<7KANJ+JI zt+^a%+&^*r{0nn@LGFViXEq?oR{(|9)&)}Ask+Stf$}~arEaN7(8})K6+l0P?gb`8w8@m1Xt8(gs#)0_UM$5=xVKM zOj>iD%ac<8>6uCG$LyBS$#cwBa@WlB>v4kVm@&hFD(9yW>XI~Hpr+OwjaU!f4*}L* zJd93=aD24KY=Z8B8|~zOk6U*YTv@p)kP=^+>++O#M)kTet?k>Wnz^8{YId~KAZ#}p zN{ufq*LNpw@TggtOWHyk0tg+$k%Phu=3Ea5PK>D4jXoMgsG3xkg&$LOMm5nBWtkyT zEa-K_5i{l0Ddn$juIsNbJM34LOx?{|SF<5twHn_^IaB$2O}D_Uz>dcce-jW@$wMHw zC&~G{xba5z;5lU9MIcg3U{gJMWPhAzpTes(_H9qiu42q^-_OsSWAsq(^+lHdLM0}( zH*Le!Ir;_9ozNA0i(egmvhL#IaXzl2v9NpUZgQZbLH`e$nPwS2xTjYWv|jqr!<8-95& zO>$b-K>gvjy>b0UqwMUtJTtR>n7YJ*E1^MFR0l(9B{9K$YV+BJebhqof;R@wweqBw z{@&E)cq5ENK0YzJ^#Lg)e`W$d_i`qj{+=JqCNsC@=Q*Y;pO4Krd*OWLUXL~JF3X!*^So8SKPAn;6Sj&U^#zB*ivpX0;^Z}7cvfYJ6m`qdC)y}^ik zwrPQIC*pR`E-QXG$>!vQSOI$GA$8{4(|;y1;*Pi6is( zow#O%>9y+;hHQH}H2bNq$@NDvj*3r9yfZOyFd$)f@z$c2-$MgN%hKGo4{iIPdqlIn zCT&q3x1;4paou(mo!?x&=kJXPX$hyM7l@M{`Vx#5?M|b?555WUFHd=6r_`k~JH5_T z8wZ_c6tA0Bta=>dxqe5>z%E&jXJW71#hWc|-+m;VJ?z8aj?q)-%sU#7haU}TG>@+H z#VWQuN%|_5>)@{0p5o=ZexqfJ?Y(^`uN_KSF1*lpC@C%RCSKbqAi-UIO!uv~S_rSl zbH=Dc6{9Z<9LVkFyg+9cv*+E=drygZJQ3sU-Gvi!lisAUyP_3(%x9gkjqRn(JMP|g zk~Euy6AwxY!x#lUK5==evb5y9JGy~T{KU!yL56j z!$!w9Z$a1NUaz56aEU*SX2EHji9A-H+9Bj_PB;?XlbDoDeq7dYq~K{# z>b6GSv-};gNmEk_Nfi}2$JbB%{?R{k+q9S7WWmb`*<428IXf23!#KV6{S%$C*U}i- zz(gaDGKKgorql7{y5V1=-e>mrJw86{H|aGCw1vAaoq7Hc8>tp*i8>YOWq5^gXN~5^ zvb?n!x~IL)w6&(|Ze!Xx44rp9?ZQkHK6T-KFD9Y?t-0IayPspf6Qp@e3ihp|OMc@! z-@C>#+@=$k>%_O%o_xNV%XowAF@F16XuTb9oH0w@2TT-ugM}@ZqTb&*+NK|FEe!ef fx;*qh=j9*7nKYij6gi3}!z(riu~v^P&tCs8A0HWJ literal 8993 zcmc&)c|6qb_E*Z1kYo=bq(n&8;oD3|wr^Bq4I$ZghB2jy>`9Sq-%p4ZWQlAumXLJ@ zA?sMeN0td=8_WDY)c4-L-+O=e{_gAdyZ86k^Loyl=RD^*&pGdN&hwe@U-h+^5AYtK zp`l^czIefqhK3fTq1i)bpoJ(&E31MKf#z4;D;m&r&z?Pd_wJ>orQNr0-~Rpk>FDU_ z>FF657#JBD4;(na#Kd&);6Y|)W)>C}R#w(ShYqo^u^m2qn4O*d$dMx)92}gSoLpR7 z+}zweJUqO-ynK9oM~@!m=jRs?5D*j;Ja+7ukdTnDu&{`T$noRHMMXu$#KgqK#ZR0# zAt51g^5n@=r%p*qN}fJ_T1rYvT3Y(dnKLpnGP1I=a&mHK&z_Z+mp^yzoPvUaqN1Xb zl9IBrvWkj|s;Vjs22)d0Q&(62<(FU1pFaa}avjE#*=OiZp{ziw)3 zYG!6;ZfTOT+1c6K+dDWo+`4tk(b4ht?b}XH zPR`EGE-o&vuC8uwZtm{x9v&W^o}OM_Uf$l`K0ZFazP^5be*XUc0RaI>Br-5CFeoSp zg+c`f2cyyGJ9qBfy?Zw#B;?+`dl(Gn{{8!*p`l@6VGkZW2oDdBh=_=cjEstkijIzs ziHV7gjeYp=VO(6?qeqY8PoE|wB|UrgEIBzjB_$;_H8m|Q zEj>N`x8HsP0O0xa=NTCpnVFecSy?Y$yvWYZ&dJGn`SRtfSFc{be*Nano7~*oyu7^p z{QS3X-xd@U6c!d16%`d17r%S=uB4=-w6yg7`}bvKW##4N6%`egm6cUhRah(*hr?A@ zSJ%|k)YjJ4)zy9Y@S(oGzM-L^v9a;v$B#`-P0h{CEiEmrt*vcsZSC#tcs!m!Aarze zbar-jb#;CE^r^eMyQin8x3~B6=g)n8ef|CY0|NttgM&juL&L+vL?UryWMp)7bZl&F ze0==NmoF0&6O)sZQ&UqQ2u@E=fBpJ(W@ct~c6M%VZhn5AL?V6r_HAKdVR3PBX=#Z} zCND29udJ-BuCA`Ft*x)GQz(>;jg8IC&8@Ai?d@$Ul{)s;0!2e3IH-N${1yL{xuJkN zM<$Z)EQH$^g|iu7kG{+<(s1GI`SUf$B96w02z$9JM=8eSWGwh{Tz+)I(){`J^Kl<& zGcMRvYFN|e2*qjebFdc)Jjgj8&f@f-;X#Pez);pF1?KMMYrSSx#zuIz)!6OMt1zo7 zvPc&}*I1FJ#GfN;nC8wJoASfG7?w;LfL7uc@RtLJ{^_*oz*byvmP(*b-PvMiws_r( zpx?DOKBgzP*P7er3w?UDYS}l?T@_pPeb=RpD|1T9WZ7L3P_I$~PYL+2@Qp|+y|3F| zRq;9{+g|q424Oy+_x7pS^fuZ&c_r=*!NFx_)8%uFmReD?XU zz6@|FwS@DhGL`ZWdq$u6mGGY>6-!sCmW}}l38jYo5f$W3k@4Lel z&>~(ddK?PolnH7wc8;W#sYE_e6oj?Wm}sJvBObCaE2!=J7D-JtTL25~W;%J1lJUBJ z;Pt+^_Lrb)u*)PaReaU6$8~vNDHv3J$TMxM=G#BdyU=ajg&ozU#K>)jCCLqL43P+= z!HKW0<{en@q1e==OZA0pw5b7(EU;$D0c=1oON(zx)oZW$yfUp8RNcC#DI(VH(KHJx zw@6aKa0)kr>z9`6j7KU7I>V--@{-oyR5>=bD!>QQJgeV25=c-HU-&X6A@ZyH_({)>Cy^oDl`caI$NKDPf~ zo~9F~dRg-FQXTwG(Tt6;BZuf39i56RIA}Z;7*QXe_sB=`&ARUY#XJb^mOAJ?i@8=| zw&xI6Izie!k2JjD`(sA0bT!5(zG8j4c0Uv6;8JKk%*AmnUQgp-k)j&l7jWZqxx6b# zO;L61%&h4;J@2Sfm+pMh=4t%`M!x)x{E8KPnq~S0|IEe5)cFi~gidKWP z!Mx*S-+#|}|9`FhcX$7Q1o_clP-%;?uj<_VnSR;1e~31qJ~b4acXBhdl;-bWc`R|W zzoz@$Ql&hg^-x_Go+w0K>vnNZDbaM(`>U@ZIN#-69~{YLgoSKcA1*%#H0Jww8jXLr z_A^tZrtt?S^4KBdW^3XTTd1_$n8htf!C705P7b!7Iprtp9sVeExUj{&uJ&7nH0uxc z$7jjefvVF3=4y*Dz{0Uw5B5Rp?G5cH?FZV_90`(a&2Y51bO-0g^Fx z%0rD@Cr$-_C%Ey2T^AE4?_#OhUwCSrexS@goF?(xTF8E;b}VcSV6b?w15FpdIgweV zcvXC7<9g}Uj{r=&z%(lV1CjjeqP{~@5Io)Vdg5pVCHkxqlSD_ zz`bZ%+;ZvcKStB>M5-hAuDeDP99R`nj6tQ54~?GIm^kywS;9W4~nvJ8}<3Kpo0bV-+OqMnav!um)(I zR2Kx3Ypqxn@n(cl$$H@;kR<0|+2t^9fe}bhYUPrvODSUobf9qckIsA03)VxC4?7D| zG6IJte+qQcXMD`2gd*3?>}3)0E*OE}$>bXZ2x43m6}ec)%{Kb`f2QLEU`shv-XMd!bqp4g z{!n)1RD}?u>AqS}Lvl<}u@}_W<~`pRli6)ATY^JF^WH!$_P94^Sgf^wHgqsG<6$)4wJ?KDfGKQ zR_?5EL5fAgEtNyaYdc5#v40h^;Zc^x)=L$@?c`sCdEk6-w`)JW+|T+aW8||ap*ySo zOh%@9VZDNC{$K|W;+dtvKKv&R)hI%J2vlF#)IgVr+9pul{U!!MB3}xqm~Yw&Vnhq@Dhw(v`gyd?X*^ z4}nAPzFGRYKvOo7Pj4+WU9j}ifTm2)v>|Ya1Dc-hBPg^oK(@%*2bn>lX2-WyRVX<# z3B|=n{%oSM^d9p+9-|N-eytSn<_2^aL$yJ}kR3amBO}E}EJ86-;@o22>9Ux-C;|st z`7>v5=mo!`OMmf>jUL~Cd{Y$SX}O6=Pd0)CCmuc+|6P%!h!I2~!w-1;&twP)DkgVp z`l+A5Am8=gc(FY4r_{_>nB&tnM)z}RkQ<-OD_-(5)3|)Ouw@2U*mn%8)RfkTBS?Pr zx_CrHui1ggmS!2~aRs$P^ydP?>{z^9?VmnTVJp7?M-XiXyVKHkm~!JG;1WYW-Hefn z|3nZR3JGq@BUx;f$hH0iOjh@n^yeHfc_KfG=+^<|bHiSAKdTbB6~daTrBD8ADe)gF zo4;dPI$;UjBledJi`>Td|DM@`;CQO)!9S$hMtyg5kNvN0E~X~0W?>R}3XAS5S+KyHBR-cHtiU;y4rIobGfp(ki2Z#NG^6bk|F5?5xR@pldb ze9~LWWEKSBkRL+b!YH%cztlEdRKj2pS=p;T5nNWnyE-)skx%E+{o|zgA2s#ghRG7+ z)>l{{5=%8_R4}==+o?dL2==J@*z^;lbS($|4c*c|Ol5FzA=ygheYtO1so3w`odB+! zIi{-z{3+_(zIV_+SMBfgU2ms+{`wPBY%dt$A8iS>kWl*focvs1i~m~oqauIQdZy7= z_(y5r8K9%n^879zNtpb$)jF@~V9~@ec;Yeg+5)cyTYg6*;`auUyr%b~SN>6SkM}~Y zvgu)p&<2m-Kxf}#XuKr)sUyo)mVQJ0&3uOg0y-e)o7=d@nBlE+UNi}-$+rUPbiKgS zjeJ}|{kHpKPnQz%2XZ^TK1Nz=#WU=9P3cEMzOyW+da&x3n~&}iZP z!=yQ5@|wiSRbc_I+D}ft%BhP5^MtV(hAd3qh zr4iiZV-Kjjcq9asz4OB44}TatS+<6>L!r>yRoJ2*P`|vTgBQS9?%VNosbvD%S=FBl zF1>fBLgWI(S^k?*A0*(0oZlB=R1AdSer(7w3`@N2z~9B_gRdaznS{0`sbk)7V~ro28C;Z)RZ+OxGb{OrEsdrR z?U+$Zn*1!h;W?G4(cf6@WMq}Mh%q5TP>Ki{t@RQrtoTpDr90`Ga zcR62hk``b7ejnGYwgX61B%G1yq_KrO?-CCs6V`^dF78~(%9it9YEicobi5-IDM3yR z$rb_)=>9ODyTN#SPK%zeQ@*pZ(wBflMV?JJr2_PwZ7KafScN(NV z22>oLwgpbzcrY3|wJeoNU`!=M4~Q~sF^?y2h0&=^E3FDt>4BW=$e8{_8=-+#ojR#i|N1nKA;j=pnBN(6Nwr$T*h9bQ7Ey0_;{D+CfGfpJ0mQYl@%tmy znAC(=CO9eb4K^DNRRgDC(JK1o)g}imT+O8bc4RT3*YZ|y!-&}xTglLUvMEZ{eF>C`TI${?L=vdFVTG?WCX5S37Wkdw-rLfBF-?`X9 zi3S^74lo7P>gvYkpby{C?A(s3D$n`|j>NX@u{7}*`Qg5zRiB2jnT9!utI$=XKXao) zey%`K_cZX=8OU%}4BeJ2=PkwaL4)qf`l13^hEpp_hLfLkj&p;jR9YQT&*FF6&PdRgi>fban2V(6i4;g& z4_a3Vz9~tUYA5N7?J}8l8&q@2EC5y8OuVBNp;uO%aF{=koY>R$eyixVZO4{py3WQ` z16vPr-(*(XN?EZHUsV0@WL6Su&=a0IAea;tNElpTkPCGivXN|cXniWwU~}a%zGcHNNtIcS)~WPmNXe z_MHOm&m&7F85?|&HHfxG<)*2}ardG$x)&O1gWgPTD5&R$L|H2v$;NuvEHzXF+l+U9 z9k4LX2#D~gC}=jh(#sdA)-V^}Y82Yd6IrLV5wsq9x8eDaoW1!6Yqf-paz1nk1=WBk z9XMa1!*_+`bp;xfKbw$Q(T;OAii0M_k$_tq48fB;tI-POlE1VTo%*W&RpFC6F?t z`O)a8Lwk{s!FrD4mwZH8A@Yjt1Keiors)N=^2FS+NymhxGj*YA0Ue?J53`HqrtUI` z<6ym!?Yg7iu1_}XEjWETE=#r6pZAK}*d$*hFUSw~Mop(TfG|EGR zMJE>|`m7>BgX~KF52|gR0`{K?$)>3_LI#cwwclQyaEs50e1Et61mD_Q7ayNbGpp}o zk9=JiyjL{95YSBNFm6guo|h2Euk}Sd%&zP}RY#8*72-HM()9yW9t1p_P;Gq+DCBv3 z!56lfRJBiIfwis~Vb|ZJI_c(LTW~B>VDlpPk(GGXc4v_<62@vO7Tbf~29EFTCpv@- z%A9)^>*}-L|CW96w8@JCwLIC>}B}_zMROyyIS~l$;IWsL9GMJ z{_{0`E*RB`N0(ugNa7b#n&*O1yY(3LSb27o`83&5nrNEJzBT>0^xo6#>+N+1x57M_ zqnF@mo$b^sy{!@MNQ)25pEUZa!)e2qBavHavm&ryqnt)BgC21)rsq5nUxQ}{j8GG*XIY4mk3UT#A=BWg|+`M!O{3;^f zkd&RGhb7&^PAP}AQNzbmC^dtQNmS9h_}0a@WtR^lk7Q$2MrtK=WC6VUZA~ewQmvP( zM(qhP0c$J$v9Y}I8^6|;L1}9lxi>v&!rCk2IQfJ<{{0(Zp-UIkNuh}FJN&Yb0d>*|8d7=e0@Y$XzU&Z@m0@hV#pBhX~8sSVOxy7xc*DtMPl auT){!>OOnj9{R5fjkbpVg*^C;kpBY8l6ai} diff --git a/vignettes/figures/WPD_histogram.png b/vignettes/figures/WPD_histogram.png index 940d978f658709e5a7fa6f7bbaa459a7db5d3332..870aa7419bec3eccfad89b83a0a87f69255f5c29 100644 GIT binary patch literal 6445 zcmb7JcUTl#vhSH8Nro#Txj{fd1yKY^(kcjufaD+p5-&L^VThwB404qri=D z8Sa%q1_@q~92X=gQ3i0xydKb9clYl5-utJ&6RPS|)vr!fb$8SaT@4mSUPb@_ShO_N z3;+NH0svx_0bT(9e)32b032P_QoCdnFt{+nVrk^i@cT0TB{6j(S1dD0UXz7M-%wkd zsXoW#wYgd4pixlqdD9Q0ue*yXM$3%~f^~)@Kc2l)P%c(+&uHr<#<B&)0s|%cK-CSrImqn0lJKIg|LRsc2xTb;pZ-f{74ZjJK{bAs0Iil88{#u zaG8c3Ui#<{o6q%TUtJxziILXdpq;Y5cu7svI7BL@ zfG_=sl%)0O##&0UT+O((NTlc1m3PY6^&mmcV!uSWqXS<#cycoe01;Km#9~CO@a!RM z6^~|%6uL&ugn6=$q$p3(KmE`YT_`T(r7ID73$<25s+TM>M7x6=%d_G2+}%2JMHUM5 z&v`BLjq0^|^M-bJ*o>RG3u`PN4G4>nt2U*pP0%MiUygG@Z)S7v=fxWIs#nZf(u2Dh zk<6w`^EXxk{o+GiLd5Xq`K3kxMn4-x?0_V@=*SOTR$eGmeh*!cQn}eI6wj5PGJnV! zFCAWf_hn0%pFLZ;)P`?Q?1<##URG|6Awjb7V?o333Mo^A8oe*?I>gG%=F)e^F@m$M z7;Ut8=F8q@A9-k7h0VP1?KavSQy8jMrm3NT*RBygasD+Y3E87f19~iV(Hk>oy-YD! z5E?=+KcR4~;YoQBvwxoLF|~gthu?f4Lhv>!K9L6GM&5B5E-MHw_P9}&MHWtS+IzvTSrf<)X3!JJ%ZUd4 zruC8%F5$u1oR8ZqdW_@>uG0Fr_CQzj8WmM19YJGjFGujB@g41I5epS}nzMQpDiTdx zRCn>@UKYX=NH-GakoLay_*cNPGG5P^*imV9!@qVoi_R^+A;0JE5uIN&rF&{_YW6@@ zCb50;4EpqFI-hB#d^&Qv!IdeSY{Dam%$H9T6b6>Jlts1kFd`vFOfpEV*{-LV;Q`*s zyC%tRIZT5^$m%n=M4t)iD30_o&RA>Ds9D)E7qHo1)@6^<*27eU%Gbu$Gx9Pd@@#}r$Yf-; z1yjS+R}Lr*dd@GBEn_S+pr|o78pYoLG<(YRA0i2AjU0NYzfMk8K4Elq1O34&M`^O?bKq|`At0Jn#2$1z0rAaU;b-l^PR zXfQ4@vHi*^zNsELX=s~)=*TT=nA=+7dM>^>vYI}Y)^Oos3{ z-!Iivl?@AsIu&@pLIB~|G)`$4BP`7VQPG`-0eU$Y74mWBu_YF`pPyp(?37l6-yaBg zGK`A-C^R;Bh85gB1op211kR^2Nnsl!qQ;z2QuHNQBQkP_%~))sY)&# z-TdM=oX0be3=pZpCj2&8U}LUJgQ};+omfAErvN(^@0xzLBK~|HEUWspaD&|@>6wIg zD%D4Z&b`|bw@s57VLKM>VAmMwxUC@XD-&VolCBgD)d^)e`|`8ER=^$>W-IMrVkaMl zhtGK5Ny94LA1wwgKDvP496_a~l~M?8_NkN6nyAWc_~*W#`l z-K0UhVHdfI^u>v4Oc8poXt?ycy#stUlD#i}?_scP=!c|AmS?LMl)PV4Q zWtihf5Z?#em2W~RCx-tiXC=?Rjn~Si^-^kUjP_Y#+^w)c*3PeG0oS z{Tj1Cd%!ft%FPJP6QZX&8qI*^o8b=~{7tz6Lj)l=q#Ar$50P=IMt_Udff<2Q!IpHm zDuMc}v0X>8{V3d=+!95*;;J5PD)| z{{miW@WA)ug5tkzUN`!*RRQz-!ALBixFT@P(AAkVI@}}j8MfllBtGy63Ztm-?d~ah zXg)3y6R`Gr7%1Hq?qxF!y1M7AE(5t4RIGfR-LU{{74UDkDDtj-=UgsA^ZGwlW$?*c zaG2NuJfbbwl*Hwf`*;doJKrAeZWV@-$k?o?Q(ZW$*X0hMRMNr6 z08)I%BvJ!e>2W%+I(_b=d#{3TF3unfTp`@4L#W3?AM zH}8RBPo1<8bX%cAl9WpA9Dnm%ri0vm2qmlx#Ok?oVEe_`R&J%rvqQlQAi^yjAZ9`) zla=iD=3UF|S#(zsS_#`P`hI1-dT$36}$!!yr|sor8H@wx_Uz&81MCkiw{Eem7&Phz`8@EKISQqrmC6Kt9y)|i4+7!JhdJ&H7I zjs;2$#itSYL4csj3i`bxDgy>g?=KbAfhbTekk!_adc3ACW~jfJpzSyV69(YK0EJ|f z=4oV9OwoyH0Gh`Sz#W1usK-3~bwxhh>Ox>S4d_SO*Gq+Th(p4q`O2x<;qWJE^iUxk z!BEYw0>GiM)afa@vr69v)PIo*M*P!~h?WTDs)y)PL-s!K;k>wP=|O;fW%)VFYvh|~ zx&D(@<$pLr0Uzjsn&BBUgY@97VtAH{7+2?;=YnVFz&fY!`g|L~I=z}P|57TY@s4sp z<|+Uv3=DXGvu4Cww0%^z96K%mR)n#UC!-yJQhmYc@y&~0fV>c}o=@(}t&yO^NTQqI z^dDPU^hi?gL%=p(xjmesMN!4)h>%&?O5NgTFkdfz38m*W*$~3^zu=wu1FkE80ky`^ zfH0$6@!sA0ikWN?WYG~hOd2(q^S8>o-XgM+uxfj7@QB2$c#TY}G-$X^-@$fsvX3D( zAzNw%U#3l zdA*-oseS2vXfsS0^qd~u_QuVU%O+x-OSaaT=2?V25xt!m>pmAD!WQThF% z>T24euEv6&f09C#F?^7!_p7haesU&uLgFg*T4OQP?K{~N{-WZ(4V4Psy^Hi9AcL%S zzc>D7bY_WIx4Fd+_flqe1J)=CFbE~c<=Lg%-)$V5`JPW3-@{U18|U^{02yja@ek+^ zFvP=n(>q%TXrN#eZT4ZM_8ACVZvQkb1ElP$wzc$vs*L4OuW;6yPi#CR6>V2jM}9P4 zDvAn>du)rkNx3vk1GlwpyAmQ_*r>w1`*K;f1d@)K-l5yyK}9pEUmF!-qq4MM{!5kk zUOTAOpIU=pnfdf*$)Rmvm`upzp9xhRfI=iTfR_O(d?o9;n^0Y#JahmA7Jg+rR94 z9_~0sauelX`wLy>cGa zvv#}42nUtR%!eiJ!~Z*Eur~qGquBPk6FOrM+;QJg35!`aMEu;x{lKLkp;^$TI9D@| zz77z*Y(6b+#+lIbv42v2i*~udzwIjp@$z&&+B;KqciTg-9`2EU^c1H6@}K?Ln{w3f z`rwP+EUU>lO6^D0*z9X|hp+CF<@pJgMpy^y4x!j7TnPxGF0GRw)GMwGR&MA|)B78(+AakWLHwJgLlD;*s|JZA~@(rq{vhTg=1G=$!O zTHD6%%+#H_aj&C$(n9Nm*WyP=80L>lx0Hu(KW@lnnVH@_EDJr4E|2fJpC9A)>H-}_ zS9wKiWKCmnKs#dg{i6ohO#6sk%O%;S2Mf$IoXn`|*qFaqdm@8_78d$4P4p0 z9UgSFPbSUPv%i@)-o9?uZfz!F*lj(@$#cEC(~gC^sLaFmTZ)DeeaN$j+d9Wsn)y1< zZ>H|{Uw5&sDtWeS*z#G`_uGXgeStmL8{LN@n9ecJMDP~B4`8}`=UumMZFaF=Ps)@UF1?EBLxj#XZRoIBO(#XM5MtzI6zH#cub=kmddq^kBcJ@`1*lNl70 z6RDS8|7X(!X`Z6Z)Z)N(qff?_XwdZ+1&j&;=G9&$~0y=h+7dA@QTQ`p{=k1?Ym<~O2W4=B>go0WCOP}X@ER!Uq1 z;Q(5JL-Pqe1iVp(!DR*Tz-<-XUgv-PB7>Rf7lBS}Ile+ogLqL$=` zDZWlXCA09hZb=F@j za^Z)7(~HM_WG8XPktDX2GyFevF!_zcd~j%c&AaH1aTQx-8zZ>#QJp-OYJFTj|8@ec zXBKa5+tUK0hR0tDi~`K;-Cf&QW38v3DPiCKGUi#eG>4lkUN|>&&o@+e^vGh<0s%=% zh!tUV*rkx;o!!f-_;A}=$&0&JU&C9PXljU~T)yAkyuk|g~PUW#u{C&vKy?3wz` zbPf&{BVr9KAP!v^-J6i3C)=p&czyAU_s(ey&f&X=z{VA=Uy}qx8r!e<1CNe3jCqaq znfKakm1+c~ShtvOsF2Ky?ONYolJmWGH1yp)U&r!C)ne~CK@&?ccY!j~FW$zZ`mozn zr=O`Fxb6KMDM2pmXi4ApS{{v0@7&k;JO9s*)0b%MuW;{Fzxcqcz`o<0MJg@nQ|s+h z;Sk_G((ChB)5xixcBBs}`n$jo1OE5Ql5o;Px^Q)?hHKB|TU-5S=M!EE_-rq{y$NvQ z4NTreb0?5`mON%fl!|JOd?IBfw%~7>ag006gf+|FnN-OIHC&_7gf>@s zop(s9j|8&$HS0eS&Y25hOhoBHL?9=*U9U85x&2wmCr-4qa^c4HJb-W#2_o{&`1r@P z4_O6XNBKmNitKp^nv)=eTouh)dafi^uoXFyH6fg2glwG_6ouCKRW_5WgoRw8L1#!X z9VBy5DMiWO9vJY;_zW*$w-aYT#Owf0EUr5XqMwXrYLpi2{%xrTP&t%WX+5br8)6F> zl*>Lu7XoXEgD0Dc4-2rT9O)ykB?TAIeb&;9!(C@m$w3)N*CABtg%x9d{)6yS0>J&p h?*+vD6Gu10DzzBqgY>RsQ~!f#sq3m0UA}eie*j|gz6$^V literal 8379 zcmch7c|4Te`@e=lg+Yrf&5%fikYt+~TZS;Q6bjjrG-Mf!wTWznsYtSw6l2Yn>=oIU z>}1P2F_y*{WB0v>r_b{|zt8ghy}sYy>-XQe&wb8y&Us(&b6xN2o(a>}(`4JvyPtu9 zflUjkZot662xVZ{v&6y(^eiToXfrT?HuZIkG=RrFd-i}pAVx+;CMKr6d-v|!w~v{b znT3Ugm6dh>{{06I95{IJAR8MSJ3BiE2L~7o=H%q$;^N}w=00@j(BZ>}j~qF2^ypC@ z9-d>zj`8yH^6~NU^YaS`2pm6t90Gwrp-@3VK_MZb6DLjx3k#n-c~V3~L{wDt)TvWq zVq&LHpFVTu%-OSN#l^)XBqSsyCC{BZCnY5%EiEl0BXj=zd0AQ63l}cP$;n;3cu`(n zUO_=YQBhGzNl95*83u#F;cx^3p`xOqs;a7{rlzj0uA!lE>Cz=lO-&>csimc*t*w3e z@?{+z9bH{rJv}`X3Z<{Ff91**0|Ns?Lqj7YBV%LZt5>gHyLQdQ#KhFp)XdE6`t|EK zZrm_8H@C2`u(Y(idGn@~m6f%%wT+F9t*xz{o!zZlx9sihZ{NP{;NWoQ&K*Zb$Gdm$ zqS0t4CnslT=X>|=xwyEvy1Kf#xw*T$V=x#G4-ZdIPcJVoEEenS?d{{^hUkPoF-0_Uu`7baYHi%=72ZV`F1qym%287x(hz%lP>CgoK2|#Kfeeq~zq} zSFc{Ze*HQnCFRYVH>s(qZ{NOs_wHRhA9D>FMe1?d|L9>+kRX{{8#Fz`)?(;Ly;}@bK`+$jFZ$KSoDK$HvCS$Hylo zCMG8*r>3TeMB>k%Kc}asXJ%$dB+~5c?A+Yk{QUgF!osg#zsO|r;^HEOLRnf`qEe~L z%gZzxZDnO;b#--ZZEbyhePd%|b8~ZRYioOZ`~1f^G6RDER!d#Q$SYxP=z#B`?ut+3 zKy?)s!JbgUUf(?_4}wF&4<9~M5q{g@(Y^|_cxcqxDFbzMk@1p2lepQZ5Bsf)itfBM zOD)oVbU!9kzsMp;HTRyMv9@CvNQ#q7KDn&S^l9ACJr582tpb^l5gJ8oAhBOQE^$0c zSk&s81A|X$xE}}+&B<~T0lmKm$FPqH!~^~xKd2ZnuBsD?x#fuI$1oW-q$!T=p`Z3$ zaN6k-xXm-im1x5XEdp!`=o$bZ0kLW*9@7jct~08RuDld$ty( zR_c4J#&T&NrmrR^q}5F7z5B{QOW+&~id~s&U}NMyx^)YqlJcPMdhaodn2 ze$?;X%pOpKh$5sk#B8L?L6zu8jhJbjLM&a^?tTYqx9KW=L@hiafZYBZ$Q&z+s%N%a z_jUQGkZRTb`p$QsFw+?q3sFpRPi*Eq=5G zg>tiu!^FM(Ww}Vn=q;+{+Cgf6O>D(@JfqmOzHy4kPxHX;(ZFRv{M~%Y>$*(M-gM^R zq|~;nl4*=3*D&2JxD9>dqYGi07xA|UtXLI+)#sIz3{S)4j@ZoO*Qu%){BU26VfOKu zmtYJ9&Knz*NN}VXL!9*Iwon|#mhF!^x8D}?<6p$38D4*$LGg&e7zR{nE4TB9dTYT` zB0H zCgRM6I`xsQ%sIxz8G9LNtDbI9=ASjiAr(tq6L%W9cHFQQ=uyp$*)Z6=>Zd~L8jDew zSOXT9SvmN$MOo0hY4e0)2D@{r$bLqd-slw zsu}z3M0{Mv*f~Gzho#*D^^6(4qG7oqtiq3e`3hIcI|sxJ8~0dbcgEq)_U8_fG`wmZ z0Mn+%YhMr0MrdrCdu3gVV3!48rEb<9ON5!wXWxr->Z6SMOA@~mXaypv)cem_S2GQ4 zsIuDPp3DXVrc_vo0M3>>DXNAhI9$Yc@n6XiB%RH|QRMvkHlE==3zx2NfVqOPq2P8&28B)YNq14Vbcm0zU zY4a?jckhQb2-JUUmo17YI+G!JV7CLQ;WzV+koGH%fiqeWzo*6dM6vs;B?3eHGevxi z^4c{N&-L!jZgA_qP#%s<6WvEd-Q==Vt*smP>&Yu_cyUaA+8%YodBi2zpKPM(%S0#k zU9+JXo%VXZ_MOpE;@5{0k{5qi3Zo8&f>5|7M*as=Qb%Dn?0)=T|)2Ycjb0 zlJdPbJtS}Yd9N{1A$ukn{{lgy5=QxKABRE&70vRT*PT3XncBJ8yMv6+;>0b<>ct7%Poj%8H`;c9eST9s#^U{)qBxjUaSK6rvdnLuHDRN`oC zmZKvx4hhA7ynV-)p}g?vAzl<9vsyg!p$!fOqKdwcOuJ zh=0o{vx{5&pG4iAZk^^pE}wRbJ#U8=rDLplmJ%F{7z?%=*5^bHzHjC`z3V4Y&D5pj zt~Kb6Z`hj3Ukl+sTE)p3F>xsZGr|zo1dr>>7j||?Sq`XNUv8LNKDK2s1vmtphlu1G zbDtnROinD)+p-wnYj^=&N-n2wK+qbO@T1goa`8)Bs~Pj`w-9viT4(qlCz~3qTRMqy z*t}s(Hv{;ULfZ;#HDt_Ww}>gBdk;A~Yj(kNZSxUt!3S=D7J2&sI)EMN4EaYh2+8a^ z^@CH7cEmr&hY{G#G-BV@AlAA44Attm#A zeeXj;xl~zfN%wm1a)!eH5s2ys>3b4(s*T7YptH=c_?iBJgW;9rQiqwiT?_^fn&3-K z70$Ws-UZYCHv!b3JmVqDmfh(dFORgv!#uSNdpvePtdVxU4Hz|uF%5zAN+CwDVbR5UbI$^um_|hQT^ceqbEaNl$2!|2g__Yryd(jO3+2 zF<{UkfDDZSp=nilIxNB_02zQ-nzU6bEr|btCNvFSf^sQkSz}lEnt~7@*$fruiz-EM zPdEBWK1n+drYrbO)p%aeJ@$-#8h~$dVTeXP?&RLG9X3BaOW8vEH>Nra!e&GMOp@ zsBtr{m+opy&*I-a2^z$eSHZ-A8R`dc$a;bLtd6IR(E$Cp#F>OXSI3$>S}&?dM5H`K z2z#-lk~7iHqWIgu?e=}<#6}0RhaO(6|K0?@K_n^K=XC;<9@GOEehIg+(Wj-S3=`kz z%2tPRSs$veeq*kP5(3NpI}`uf@IUu*Y08h}K14iFg}RJSZ=1${KIC+Hc{Pk6xH?IfX-hANTe0 zt~vr@)xHeH9LQBlU--6&9p*R**78ZB$us+tdd=v3bcM z`Nu2ic>M6% zIG$tGV%R6z8ygCo6ejGk0@_g5gOtzG&gmZ8{I+vwZYQjGl^Ux$-U zPRklDFX^0Wp~qNN1^gp>blW<8KM7F}H}Sn8|=4fTBW&JAZP+Z>{3|##}l3 zkQ0LSU#Cm&VCXs?@tT_GBf+c50N2MfGc7Lz25@|7i)sb z-F5!jm!C!IKO76$DEo7AzRzJ&AhUDT!q5g_Pz6ReD6}|G&SY>9ps_&xw_7{#ZnQE3 zu|S}{G3Jb=|IaSx3Ct;6Jy{U~q;z7e37Z{#&==okV9?iN^<_i!^9QZI0Frw6xDUQ< z0e&(U@JHmxDwOU`f|gn0$ZIB{x&>b#@a$ewBy-1`6w_i_IRtR#@70RiuCGSHqe@9{ zqCC9e0HiC4fXT-BHYh}ND;=)%r_ot0?)~R>^6ER z)8tGpgc6~#3helSFa1~=6B>LhHE@eUL0+%+VC0wvU&@1{=ulyx7WD>D#!1@SD>I}M@V4eh-hF!JV2pGvK$gX8q;B~dD5R}-;&e?Wvjak~9*7g&BQ zw%A!GPPu@p0{G72=oxxWIE<@dZz{JE!h78r_bl2~5K~RTE!ZN@b;E)!Wr-cZCBJZ) z5&~Jw^vvPg$;35scae(u&Wv zdt-8UGJhp5I@@EojD(}V+L-v!5C$jmD=Y7&u|Cs+bL6_mY5D&)um2^*|EK1!BgodH zr`=BAt$|QY&snN=PP>E#vXn)7a|9wBhY+4%ujyRB@{toM1r7yrf|oaZzZ^(xALt2Q z5sZ>Fnz{u23Md62kJ(k|_wOW~ytET6sXQ#|59*Qw?-JYwR|httwK|uVZvQb2e}&Dv znl)pyBwohz*U77eJ9ht;>~@mM{?u|Qaw4+$cgj(N@=-rUnB|fZfv`g#F~L4jIRAF? z1a;A=0MKwqzz~yUrU@JA>U8W;TeQedI@%~L19a>IaonH_>5QK9;jKdJM+o6QE+q9E zSugF#R^b?wSk)$of7QoTuTfm)!G31LK_f9OAXyF2s`!uOr3E)2cK`xAY0Hlt$Mg^} zR@=EV9*i9HTBFJDsP)9GfaGUDHfIrz1bD)RtNEJ>lc1680A*x*U#(oO64u zmb)#nA+K5U!5n>87|iOU0j)Tqu&_B3{3Y3qsS8EwREb$iDn0Rs^gAY3`tC!XhX&Q| z5ACnJQ?c6bTtt1SEG4mW{v{a1+Vy8i5Og0T@9?zqNe-(Z9+5rUUgiCzFJnIU7QK?# zv)7tF`fBGGwv#B~(D+^NpPkxvi95sI;{e@gvyu+G*L4&_kTQx1~ zq)zESN|WC;CB15C9rFjY7QG-*jIc!eDE_mq|8LOuVz;%;?(7*$bm0ezyz7thm68EH`5ba`<&V%Y=OQmfA~%9ZK0uS8D1N3_}S zuh*di9z}Ui5aI8T6qm`xZ%J0R6YQz}X36XsskAPEIK|}$h&!bvTcA0`K@82Ican#= z)wHa;?Z#cQ`f-K5cYiIP0ZK~Sm08-!T^DUL7e3P@?K!&_ch(J_+X$Cm7~1yhoxB%N zpSd*=;FKZ<8BaLF(|a&++&z;tJXm#Ec^f~;o|*!~baoZkHIZ6+e=GvZj(4owwOBy5 z2kYE%B`49&Fo6~vac^#xRNC7qD`=VPIybwbOsx1!Eb}-%TI%_&yVd6n$|GvL)U(j> zu%j!aw>6(HfmZ9=PED?9&k12K{aoaxaZlz*;_4}QEun%mjI)zazCkP9_=H8zU;I>5 zP_X{`St3Wp_N89qifI_J;usJV^<0MdNU-F1CCm9i+-^7bli;`d1bp0^;ng+M+!$y(1 z0=BrT)rQiJPh30EtV!yBcn(s+_>U@Wd`U(oUJCMh@qKuAq1W(I$>o<5ZF`f>L97s#yu_g`sxj_><{c-&iW{j;nv*#kUg|RT4VBG^R0@sVA*$-i$o( zv$DaBzmA(G)A}?+&4LPV&2dZWAI!T?Jgg3Ns}oT;>p)5jqtGJ1k9mepZ0vKjG8i<{R+nVOJ zorNFSmyd}7!ttRxwYu9C?>ld22e^mlcxePI{)9PVR(UHnW6dqkEE6WbCLuDCI^lx> zX|TDxYx_vDX*W91^&K@6?(kV=;4`gr(Yn4RZ7wIGbWO1VQ3%^0+d|Do49sT@n;t3S zbSgChN@m*>714V8Qf|$c+N4`tE%8|I2R_AyNOWlC%2!^3tcBP@6^X^`m))R}!Bo(* z&+kR?_NZRX8;Wkm0f`4~tLD7XaKAaUSuYsra$~Kc>e`Ih$beD4`L)^Jyon$E0DM%H z@&*mBS=Q8?yVVT`Bua^EMo`U@A}PtoDerq0xREYO?mgAOFFhd}Q^G?VC!gyv=O?s5 zHKg4;QUfp$iq2_#B>|CBW=