Wind-Energy-Yield-Forecast.qmd

---
title: "Wind Power & Energy Yield Forecast"
subtitle: Forecast Horizon -> Seven days ahead
author: "Jorge A. Thomas"
date: "`r Sys.Date()`"
format:    
    html:
      self-contained: true
      code-fold: true
      df-print: tibble
      code-summary: "Show the code"
      grid: 
        margin-width: 350px
execute: 
  echo: fenced
reference-location: margin
citation-location: document

---

Wind power and on-site energy forecast for **a Goldwind turbine model GW 165-6.0 6300**.

## Location and Site Assessment

```{r}
#| label: init
#| message: false
#| echo: false

# ======================== Wind-Energy-Yield-Forecast =========================
# June 2024
# @jorgethomasm

# ========================== Load Dependencies ==========================

library("reticulate")

# Use a specific virtualenv
# reticulate::use_virtualenv(Sys.getenv("VIRTUAL_ENV"), required = TRUE)
reticulate::use_virtualenv(paste0("/opt/", Sys.getenv("VENV_NAME")), required = TRUE)

# Sys.getenv("VIRTUAL_ENV")
# Sys.setenv(LD_LIBRARY_PATH = paste0(Sys.getenv("HOME"),"/opt/python-3.6.2/lib"))
# Sys.getenv("LD_LIBRARY_PATH"
# py_config()

library("dplyr")
library("ggplot2")
library("tibble")
library("lubridate")
library("kableExtra")
library("xts")
library("plotly")

theme_set(theme_minimal())

source("./utils/jthomfuncs_wind.r")
```

::: {#fig-intro layout-ncol="1"}
![Source: Google Earth](./imgs/wt_location.png){width="100%"}

Location of the wind turbine.
:::

```{r}
#| label: calculations
#| message: false
#| echo: false

reticulate::source_python("./utils/request_weather.py")

hourly_df <- py$hourly_dataframe

# Last update weather data logSys
write.table(hourly_df, "./data/last_req_weather.csv", sep = ",", quote = FALSE)

# ====== Wind Turbine Specs ======

# Manufacturer: Goldwin
# Model: GW 165-6.0 6300 

# Variable Speed - Variable pitch

wt_rated_power <- 6000 #  [kW]
wt_rated_wind_speed <- 11.4 # [m/s] at standard air density
wt_area <- 21382 # squared metres
wt_hub_height <- 100 # metres
Cp <- 0.45 # Power Coefficient
wt_rotor_diam <- 165 # metres

wt_ws_cut_in <- 3 # m/s
wt_ws_cut_out <- 24 # m/s

# Speed Range During Power Production
wt_min_speed <- 5.5  # [RPM]
wt_max_speed <- 10.7  # [RPM] Nominal

#' Tip speed of blade 
#' Linear speed of blade tip for Tip-Speed Ratio (lambda) calc.
wt_min_tip_speed <- 2*pi * (wt_min_speed/60) * (wt_rotor_diam/2) # [m/s]
wt_max_tip_speed <- 2*pi * (wt_max_speed/60) * (wt_rotor_diam/2) # [m/s] Nominal

# wt_max_tip_speed <- 92 # [m/s] # From Specs.

# ====== Onshore Wind Park =======

n_wind_turbines <- 1
installed_capacity <- (n_wind_turbines * wt_rated_power) / 1000 # [MW]


 # Lat =  53°52'44.10"N
 # Lon =   7°24'28.35"E


# ====== Power & Energy Calculations =======

# Add calculations:

hourly_df <-
	hourly_df |>
  mutate(date =  lubridate::with_tz(date, "UTC")) |> 
	mutate(rho_h = calc_humid_air_density(temperature_120m, relative_humidity_2m/100, surface_pressure)) |>
  mutate(p_out = calc_wt_output_power(rated_power = wt_rated_power,
                                        area = wt_area,
                                        power_coeff = Cp,
                                        cut_in = wt_ws_cut_in,
                                        cut_out = wt_ws_cut_out,
                                        air_density = rho_h,
                                        wind_speed = wind_speed_120m)) 

# To XTS

hourly_xts <- xts::xts(hourly_df[, 6:ncol(hourly_df)], as.POSIXct(hourly_df$date, tz = "UTC"))

ends <- xts::endpoints(hourly_xts, "days")
#ends <- ends + 1
#ends[length(ends)] <- ends[length(ends)] - 1

```


## Forecasted Generated Power by the Wind Turbine
```{r}
#| label: plots
#| echo: false
#| fig-width: 11
#| fig-height: 5
#| column: page-right
#| warning: false

# library(extrafont)
# loadfonts(quiet = TRUE)

pl_power <- ggplot(hourly_df, aes(x=date, y=p_out)) +

  geom_line(colour = "red") + geom_point(colour = "red") +
  labs(title = strftime(Sys.Date(), format = "%Y", tz = "UTC") ,y="[kW]", x="UTC", caption = "Forecast Horizon: 7 days ahead.") +
  scale_x_datetime(date_breaks = "1 day", date_minor_breaks = "2 hour", expand = c(0, 0), date_labels = "%d.%b")

ggplotly(pl_power)

```

## Daily Forecasted Energy Yield of the Wind Turbine

```{r}
#| label: daily_energy
#| echo: false
#| message: false
#| fig-width: 11
#| fig-height: 5
#| column: page-right
#| warning: false

daily_xts <-  xts::period.apply(hourly_xts, ends, sum)

daily_df <- 
    daily_xts |>
    as_tibble() |>
    mutate(datetime = as.Date(time(daily_xts), tz = "UTC")) |>
    mutate(MWh = V1/1000) |>
    relocate(datetime) 
   
p_eng <- ggplot(daily_df, aes(x=datetime, y=MWh)) +
  geom_col(alpha = 0.6, fill="blue") + 
  labs(y = "[MWh]", x="Days Ahead", caption = "Forecast Horizon: 7 days ahead.") +
  geom_text(aes(label = round(MWh, 2)), vjust = -0.2) +
  scale_y_continuous(limits=c(0, 200)) +
  scale_x_date(date_breaks = "1 day", date_minor_breaks = "1 day", expand = c(0, 0), date_labels = "%d.%m.%Y")

ggplotly(p_eng)

```

#### Requested Weather Variables (DWD ICON)

- relative_humidity_2m \[%\]
- surface_pressure \[hPa\] 
- wind_speed_120m \[m/s\] 
- temperature_120m \[°C\]

#### Tech Stack

- CI/CD: Docker & GitHub Actions
- Quarto
- Python
- Open Meteo
- Numeric Weather Prediction (NWP): DWD ICON
- R
- ggplot
- plotly

::: {#fig-intro layout-ncol="1"}

![GW 165-6.0 6300](./imgs/Goldwind_GW165-60.jpg){width="35%"}

Wind turbine model.
:::