Code
library(shiny)
library(dplyr)
library(lubridate)
library(ggplot2)
library(glue)solar wind model
Set up
Plot function
Code
plot_vre_correlation <- function(cf_solar,
cf_wind,
phi_target,
solar_share,
wind_share,
scenario = "",
dt_min = 1,
days = 5) {
# 1. Day/night
sunrise <- 9
sunset <- 17
day_len <- sunset - sunrise
night_len <- 24 - day_len
# Clamp CF to theoretical max
cf_solar <- min(cf_solar, day_len/24)
cf_wind <- min(cf_wind, night_len/24)
# 2. Raw ramps → clamp → plateau
raw_ramp_s <- day_len - 24*cf_solar
ramp_s <- min(raw_ramp_s, day_len/2)
plat_s <- day_len - 2*ramp_s
raw_ramp_w <- night_len - 24*cf_wind
ramp_w <- min(raw_ramp_w, night_len/2)
plat_w <- night_len - 2*ramp_w
# 3. Time setup
dt_h <- dt_min/60
steps_per_day <- 24*60/dt_min
n_steps <- steps_per_day * days
times <- seq(
from = as.POSIXct("2025-01-01 00:00"),
by = paste0(dt_min, " min"),
length.out = n_steps
)
eff <- dt_h/2
sr_eff <- sunrise - eff
ss_eff <- sunset + eff
# 4. Build raw trapezoids
df <- tibble(datetime = times) %>%
mutate(
h24 = (hour(datetime) + minute(datetime)/60) %% 24,
solar_raw = case_when(
h24 < sr_eff ~ 0,
h24 < sr_eff + ramp_s ~ (h24 - sr_eff) / ramp_s,
h24 < ss_eff - ramp_s ~ 1,
h24 < ss_eff ~ (ss_eff - h24) / ramp_s,
TRUE ~ 0
),
h_rel = (h24 - sunset + 24) %% 24,
wind_raw = pmax(0, pmin(
(h_rel - eff) / ramp_w,
(night_len + eff - h_rel) / ramp_w,
1
))
)
# 5. Normalise to CF×24
A_s <- sum(df$solar_raw[1:steps_per_day]) * dt_h
A_w <- sum(df$wind_raw[1:steps_per_day]) * dt_h
df <- df %>%
mutate(
solar = solar_raw * ((24*cf_solar)/A_s),
wind = wind_raw * ((24*cf_wind)/A_w)
) %>%
# Apply your 2-unit capacity mix
mutate(
solar = solar * (2 * solar_share),
wind = wind * (2 * wind_share)
)
# 6. Continuous Pearson lag (find the best lag)
sun <- df$solar[1:steps_per_day]
win <- df$wind[1:steps_per_day]
lags <- 0:(steps_per_day-1)
phis <- sapply(lags, function(L) {
wsh <- if (L == 0) win else c(win[(L+1):steps_per_day], win[1:L])
cor(sun, wsh, use = "pairwise.complete.obs")
})
best <- lags[which.min(abs(phis - phi_target))]
# 7. Shift wind
wsh <- if (best > 0) c(win[(best+1):steps_per_day], win[1:best]) else win
df$wind <- rep(wsh, days)
# 8. Coverage metrics
combo <- df$solar + df$wind
demand <- max(combo)
met <- sum(pmin(combo, demand)) * dt_h
pct_m <- met / (demand * days * 24) * 100
caption_text <- if (nzchar(scenario)) {
scenario
} else {
""
}
subtitle <- sprintf("Solar + wind = %.1f%% demand", pct_m)
# 9. Plot
df %>%
mutate(hour = as.numeric(difftime(datetime, first(datetime), units = "hours"))) %>%
ggplot(aes(x = hour)) +
geom_area(aes(y = wind), fill = "darkgreen", alpha = 0.8) +
geom_ribbon(aes(ymin = wind, ymax = wind + solar), fill = "gold", alpha = 0.8) +
geom_hline(yintercept = demand, linetype = "dashed", color = "red") +
scale_x_continuous(breaks = seq(0, 24*days, 4),
labels = function(x) sprintf("%d", x %% 24)) +
scale_y_continuous(limits = c(0, NA), expand = c(0,0)) +
labs(
title = "Variable Renewable Energy Penetration",
subtitle = subtitle,
caption = caption_text,
x = "Hour of Day",
y = "Power (units)"
) +
theme_minimal() +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size = 14),
plot.subtitle = element_text(size = 12)
)
}Manual print
Code
plot_vre_correlation(0.25, 0.35, -1, 0.5, 0.5)
Code
plot_vre_correlation(0.25, 0.35, 1, 0.5, 0.5)
Code
plot_vre_correlation(0.25, 0.35, -0.2, 0.5, 0.5, "-0.2 wind and solar correlation")
Code
plot_vre_correlation(0.25, 0.35, -0.2, solar_share = 0.74, 0.26, "-0.2 correlation and 2.8MW of solar capacity for every MW of wind")
App
Code
# — UI & Server —
ui <- fluidPage(
titlePanel("Interactive VRE Coverage Simulator"),
sidebarLayout(
sidebarPanel(
width = 3,
# Solar CF
fluidRow(column(12,
numericInput("cf_solar_num", "Solar CF", 0.25, min = 0, max = 1, step = 0.01),
sliderInput("cf_solar", NULL, 0, 1, 0.25, step = 0.01)
)),
# Wind CF
fluidRow(column(12,
numericInput("cf_wind_num", "Wind CF", 0.35, min = 0, max = 1, step = 0.01),
sliderInput("cf_wind", NULL, 0, 1, 0.35, step = 0.01)
)),
# Correlation
fluidRow(column(12,
numericInput("rho_num", "Correlation (ρ)", -1, min = -1, max = 1, step = 0.01),
sliderInput("rho", NULL, -1, 1, -1, step = 0.01)
)),
# Solar share
fluidRow(column(12,
numericInput("s_share_num", "Solar share", 0.5, min = 0, max = 1, step = 0.05),
sliderInput("s_share", NULL, 0, 1, 0.5, step = 0.05)
)),
# Wind share
fluidRow(column(12,
numericInput("w_share_num", "Wind share", 0.5, min = 0, max = 1, step = 0.05),
sliderInput("w_share", NULL, 0, 1, 0.5, step = 0.05)
))
),
mainPanel(
plotOutput("vrePlot", height = "500px"),
tags$div(style = "text-align: center; margin-top: 10px;",
actionButton("refresh", "Refresh"),
actionButton("reset", "Reset")
)
)
)
)
server <- function(input, output, session) {
# Sync sliders ↔ numeric inputs
observeEvent(input$cf_solar, updateNumericInput(session, "cf_solar_num", value = input$cf_solar))
observeEvent(input$cf_solar_num, updateSliderInput( session, "cf_solar", value = input$cf_solar_num))
observeEvent(input$cf_wind, updateNumericInput(session, "cf_wind_num", value = input$cf_wind))
observeEvent(input$cf_wind_num, updateSliderInput( session, "cf_wind", value = input$cf_wind_num))
observeEvent(input$rho, updateNumericInput(session, "rho_num", value = input$rho))
observeEvent(input$rho_num, updateSliderInput( session, "rho", value = input$rho_num))
observeEvent(input$s_share, updateNumericInput(session, "s_share_num", value = input$s_share))
observeEvent(input$s_share_num, updateSliderInput( session, "s_share", value = input$s_share_num))
observeEvent(input$w_share, updateNumericInput(session, "w_share_num", value = input$w_share))
observeEvent(input$w_share_num, updateSliderInput( session, "w_share", value = input$w_share_num))
# Reset button: revert all controls to defaults
observeEvent(input$reset, {
updateNumericInput(session, "cf_solar_num", value = 0.25)
updateSliderInput( session, "cf_solar", value = 0.25)
updateNumericInput(session, "cf_wind_num", value = 0.35)
updateSliderInput( session, "cf_wind", value = 0.35)
updateNumericInput(session, "rho_num", value = -1)
updateSliderInput( session, "rho", value = -1)
updateNumericInput(session, "s_share_num", value = 0.5)
updateSliderInput( session, "s_share", value = 0.5)
updateNumericInput(session, "w_share_num", value = 0.5)
updateSliderInput( session, "w_share", value = 0.5)
})
# Only re-run when Refresh is clicked (and once at launch)
params <- eventReactive(input$refresh, {
list(
cf_solar = input$cf_solar,
cf_wind = input$cf_wind,
phi_target = input$rho,
solar_share = input$s_share,
wind_share = input$w_share
)
}, ignoreNULL = FALSE)
output$vrePlot <- renderPlot({
p <- params()
plot_vre_correlation(
cf_solar = p$cf_solar,
cf_wind = p$cf_wind,
phi_target = p$phi_target,
solar_share = p$solar_share,
wind_share = p$wind_share,
dt_min = 1,
days = 5
)
})
}
shinyApp(ui, server)Shiny applications not supported in static R Markdown documents