Monte Carlo Simulation and Bootstrap for Estimating the Cost of Wind Energy

Monte Carlo Simulation and Bootstrap for Estimating the Cost of Wind Energy

Objective Description

This analysis estimates the Levelized Cost of Wind Energy (LCOE) using Monte Carlo simulation to account for uncertainties in key parameters (investment cost, operation and maintenance cost, interest rate, capacity factor, and lifetime) and bootstrap methods to compute confidence intervals for the mean LCOE. It includes wind speed modeling, vertical extrapolation, power density calculation, Weibull distribution fitting, and correlation analysis.

Packages

To install these packages, use: install.packages('package_name').

library(psych)
library(Hmisc)
library(MASS)
library(fitdistrplus)
library(corrplot)
library(scatterplot3d)
library(plot3D)

Wind Speed Data Generation

Generate synthetic wind speed data with normal distribution for eastward and northward components.

set.seed(42)
n <- 100
u10 <- rnorm(n, mean = 0, sd = 8)   # Eastward wind (m/s)
v10 <- rnorm(n, mean = 0, sd = 8)   # Northward wind (m/s)
data <- data.frame(u10, v10)
windSpd <- function(u10, v10) {
  sqrt(u10^2 + v10^2)
}
data$wind_speed <- windSpd(data$u10, data$v10)
if (any(is.na(data$wind_speed)) || any(data$wind_speed <= 0)) {
  stop("Invalid wind speed data: contains NA or non-positive values")
}
head(data)

##          u10        v10 wind_speed
## 1 10.9676676  9.6077230  14.580743
## 2 -4.5175854  8.3580087   9.500783
## 3  2.9050273 -8.0256692   8.535253
## 4  5.0629008 14.7878552  15.630535
## 5  3.2341466 -5.3341873   6.238049
## 6 -0.8489961  0.8441105   1.197212

Vertical Extrapolation of Wind Speed

Extrapolate wind speed from 10 meters to 80 meters using a power-law model.

v1 <- data$wind_speed
h1 <- 10   # Base height (m)
h2 <- 80   # Target height (m)
alpha <- (0.37 - 0.088 * log(mean(v1))) / (1 - 0.088 * log(h1/10))
v2 <- v1 * ((h2/h1)^alpha)
if (any(is.na(v2)) || any(v2 <= 0)) stop("Invalid extrapolated wind speed")
mean(v2)

## [1] 13.74366

Wind Power Density

Calculate wind power density using air density and cubed wind speed.

rho <- 1.225 # kg/m^3
PDA <- (rho * mean(v2^3)) / 2
PDA

## [1] 3280.681

Weibull Distribution Fitting

Fit a Weibull distribution to the extrapolated wind speed using multiple methods.

fit.weibull <- fitdist(v2, distr = "weibull", method = "mle", lower = c(0, 0))
fit.weibull2 <- fitdist(v2, distr = "weibull", method = "mse", lower = c(0, 0))
fit.weibull3 <- fitdist(v2, distr = "weibull", method = "mge", lower = c(0, 0))

## Warning in fitdist(v2, distr = "weibull", method = "mge", lower = c(0, 0)):
## maximum GOF estimation has a default 'gof' argument set to 'CvM'

gofstat(list(fit.weibull, fit.weibull2))

## Goodness-of-fit statistics
##                              1-mle-weibull 2-mse-weibull
## Kolmogorov-Smirnov statistic    0.07193344    0.08373539
## Cramer-von Mises statistic      0.04310802    0.07047810
## Anderson-Darling statistic      0.34143359    0.49751189
## 
## Goodness-of-fit criteria
##                                1-mle-weibull 2-mse-weibull
## Akaike's Information Criterion      679.4113      679.8889
## Bayesian Information Criterion      684.6217      685.0993

qmedist(v2, "weibull", probs = c(1/3, 2/3))$estimate

##     shape     scale 
##  1.755781 15.633745

# Extract Weibull parameters for verification
k <- fit.weibull$estimate["shape"]
c_scale <- fit.weibull$estimate["scale"]
print(c(shape = k, scale = c_scale))

## shape.shape scale.scale 
##    1.876741   15.459945

Capacity Factor Calculation

Calculate the capacity factor based on Weibull parameters and turbine characteristics.

k <- 1.755781    # Shape parameter
c_scale <- 15.633745   # Scale parameter
vi <- 4; vr <- 10; vo <- 25  # Turbine characteristics (m/s)
Pout <- exp(-(vi/c_scale)^k) - exp(-(vr/c_scale)^k)
Pr <- (vr/c_scale)^k - (vi/c_scale)^k
if (Pr == 0) stop("Division by zero in capacity factor calculation")
CF <- Pout/Pr - exp(-(vo/c_scale)^k)
CFW <- (exp(-(vi/c_scale)^k) - exp(-(vr/c_scale)^k)) / ((vr/c_scale)^k - (vi/c_scale)^k) - exp(-(vo/c_scale)^k)
NominalPower <- 5200 # kW
TimeW <- 20  # Years
Egw <- 8760 * TimeW * NominalPower * CFW
round(Egw, 3)

## [1] 603498909

LCOE Estimation

Estimate the Levelized Cost of Energy (LCOE) for the wind turbine.

UPw <- 1600        # $/kW
instCost <- 0.3    # Installation cost (30%)
COMw <- 0.25       # O&M cost (25%)
it <- 0.087        # Interest rate
NominalPower <- 5200  # Ensure defined
TimeW <- 20  # Ensure defined
if (!exists("Egw")) stop("Egw not defined. Ensure capacity_factor chunk runs successfully.")
LCOEwind <- ((NominalPower * UPw + instCost * (UPw * NominalPower)) * 
             (1 + COMw * (((1 + it)^TimeW - 1) / (it * (1 + it)^TimeW)))) / Egw
round(LCOEwind, 3)

## [1] 0.06

Monte Carlo Simulation

Simulate uncertainties in key parameters to estimate LCOE distribution.

set.seed(42)
n <- 1000  # Reduced for testing
IC_mean <- 2496000
IC_uncert <- 0.5
IC_sd <- IC_mean * IC_uncert
IC_normal <- rnorm(n, mean = IC_mean, sd = IC_sd)
hist(IC_normal, prob = TRUE, breaks = 50, col = "blue",
     xlab = "Investment Cost ($/kW)", ylab = "Density",
     main = "Monte Carlo Simulation - Investment Cost")
box()

OM_mean <- 0.25
OM_uncert <- 0.5
OM_sd <- OM_mean * OM_uncert
sOM <- rnorm(n, mean = OM_mean, sd = OM_sd)
hist(sOM, prob = TRUE, breaks = 50, col = "green",
     xlab = "O&M Cost ($)", ylab = "Density",
     main = "Monte Carlo Simulation - O&M Cost")
box()

IR_mean <- 0.1087
IR_uncert <- 0.5
sIR <- runif(n, min = IR_mean * (1 - IR_uncert), max = IR_mean * (1 + IR_uncert))
hist(sIR, prob = TRUE, breaks = 50,
     xlab = "Interest Rate (%)", ylab = "Density",
     main = "Monte Carlo Simulation - Interest Rate")
box()

CF_mean <- CFW
CF_uncert <- 0.3
sCF <- runif(n, min = CF_mean * (1 - CF_uncert), max = CF_mean * (1 + CF_uncert))
sCF[sCF > 1] <- 1
sCF[sCF < 0] <- 0
hist(sCF, prob = TRUE, breaks = 50, col = "blue",
     xlab = "Capacity Factor", ylab = "Density",
     main = "Monte Carlo Simulation - Capacity Factor")
box()

LT_mean <- 20
LT_uncert <- 0.3
sLT <- runif(n, min = LT_mean * (1 - LT_uncert), max = LT_mean * (1 + LT_uncert))
hist(sLT, prob = TRUE, breaks = 50, col = "orange",
     xlab = "Lifetime (years)", ylab = "Density",
     main = "Monte Carlo Simulation - Lifetime")
box()

sLCOE <- (IC_normal * (1 + (sOM * (((1 + sIR)^sLT) - 1) / (sIR * (1 + sIR)^sLT)))) /
         (NominalPower * sCF * 8760 * sLT)
summary(sLCOE)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.005953  0.006704  0.011860  0.013199  0.017822  0.090258

hist(sLCOE, breaks = 50, prob = TRUE, border = "yellow", col = "blue3",
     main = "Monte Carlo Simulation - LCOE", xlab = "LCOE ($/kWh)", ylab = "Density")
lines(density(sLCOE), col = "red", lty = 1, lwd = 1)

Monte Carlo Results

Summarize and quantify the LCOE distribution.

round(cbind(quantile(sLCOE, c(0.05, 0.50, 0.95))), 5)

##        [,1]
## 5%  0.00126
## 50% 0.01186
## 95% 0.03051

library(scales)

## 
## Attaching package: 'scales'

## The following objects are masked from 'package:psych':
## 
##     alpha, rescale

cbind(scientific(quantile(sLCOE, c(0.05, 0.50, 0.95)), digits = 3))

##     [,1]      
## 5%  "1.26e-03"
## 50% "1.19e-02"
## 95% "3.05e-02"

round(cbind(median(IC_normal), median(sOM), median(sIR), median(sCF), median(sLT)), 3)

##         [,1]  [,2]  [,3]  [,4]   [,5]
## [1,] 2479609 0.249 0.111 0.666 19.981

Correlation Analysis

Analyze correlations between simulated parameters and LCOE.

data_mc <- data.frame(IC_normal, sOM, sIR, sCF, sLT, sLCOE)
mcor <- cor(data_mc)
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(mcor, method = "color", col = col(200), type = "full",
         addCoef.col = "black", order = "hclust")

Bootstrap Analysis

Estimate confidence intervals for the mean LCOE using bootstrap resampling.

set.seed(42)
n_boot <- 100
mean_LCOE <- numeric(n_boot)
for (i in 1:n_boot) {
  sample_i <- sample(sLCOE, size = length(sLCOE), replace = TRUE)
  mean_LCOE[i] <- mean(sample_i)
}
summary(mean_LCOE)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01258 0.01293 0.01313 0.01316 0.01332 0.01441

quantile(mean_LCOE, probs = c(0.025, 0.975))

##       2.5%      97.5% 
## 0.01267612 0.01384998

3D Scatterplot

Visualize relationships between investment cost, interest rate, and LCOE.

set.seed(42)
n_points <- 50
x <- sample(IC_normal, n_points, replace = FALSE)
y <- sample(sIR, n_points, replace = FALSE)
z <- sample(sLCOE, n_points, replace = FALSE)
scatter3D(x, y, z, pch = 16, cex = 1.5,
          colvar = z, col = ramp.col(c("blue", "yellow", "red")),
          xlab = "Investment Cost ($/kW)", ylab = "Interest Rate", zlab = "LCOE ($/kWh)",
          main = "3D Scatterplot Monte Carlo",
          colkey = list(length = 0.5, width = 1, cex.axis = 0.8, cex.clab = 0.8, side = 4, ticktype = "simple"))

## Warning in overrulepar(parameter, colkeypar): unknown names in colkey parameter
## subset: ticktype

## Warning in overrulepar(parameter, colkeypar): unknown names in colkey parameter
## subset: ticktype

## Warning in axis(side = 4, mgp = c(3, 1, 0), las = 2, at = NULL, labels = TRUE,
## : "ticktype" is not a graphical parameter

Results Extraction

Save LCOE results to a file.

write.table(sLCOE, "LCOE_results.txt", col.names = FALSE, row.names = TRUE)
getwd()

## [1] "C:/Abdi-Basid ADAN/22.Doc.Rpubs/2 Rpubs SM LCOE H A W"

```