This analysis uses:
See below
We save this locally, if it is older than 1 day we re-download.
# this needs to be more clever - only download dates we want from API?
esoF <- paste0(esoDataPath, "latest_df_fuel_ckan.csv")
# check for gzipped version (see below)
if(file.exists(paste0(esoF, ".gz"))){
message("We already have a version saved to: ", paste0(esoF, ".gz"))
message("Loading it...")
ngeso_dt_raw <- data.table::fread(paste0(esoF, ".gz"))
} else {
message("We don't already have a version, downloading and saving to: ", esoF)
ngeso_dt_raw <- data.table::fread("https://data.nationalgrideso.com/backend/dataset/88313ae5-94e4-4ddc-a790-593554d8c6b9/resource/f93d1835-75bc-43e5-84ad-12472b180a98/download/df_fuel_ckan.csv")
data.table::fwrite(ngeso_dt_raw, esoF) # save locally for future re-use
dkUtils::gzipIt(esoF)
}
## We already have a version saved to: /Users/ben/Dropbox/data/UK_NGESO/genMix/latest_df_fuel_ckan.csv.gz
## Loading it...
# if older than 1 day, reload
today <- lubridate::today()
lastNGESO <- as.Date(max(ngeso_dt_raw$dv_start))
if(today - lastNGESO > 1) {
# old data, reload
message("But the version we have dates from ", lastNGESO, " (",today - lastNGESO ," days ago), downloading latest...")
ngeso_dt_raw <- data.table::fread("https://data.nationalgrideso.com/backend/dataset/88313ae5-94e4-4ddc-a790-593554d8c6b9/resource/f93d1835-75bc-43e5-84ad-12472b180a98/download/df_fuel_ckan.csv")
# nice dateTime
ngeso_dt_raw[, dv_start := lubridate::as_datetime(DATETIME)]
data.table::fwrite(ngeso_dt_raw, esoF)
dkUtils::gzipIt(esoF)
}
# we think renewable is wind + solar, low carbon includes nuclear
ngeso_dt <- addDerivedVariables(ngeso_dt_raw, source = "NGESO")
ngEsoCap <- paste0("Source: NG-ESO generation mix data ",
min(ngeso_dt$dv_start), " - ",
max(ngeso_dt$dv_start))
# check
#table(ngeso_dt$dv_peakPeriod)
ggplot2::ggplot(ngeso_dt, aes(x = dv_month,
group = dv_month,
y = CARBON_INTENSITY)) +
geom_boxplot() +
labs(x = "Month",
y = "Halfhourly carbon intensity (g CO2/MW)",
caption = ngEsoCap)
plotDT <- ngeso_dt[, .(meanGen_GW = mean(GENERATION/1000)),
keyby = .(dv_hh_start_hms, dv_year)]
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_hms,
y = dv_year,
fill = meanGen_GW)) +
geom_line() +
labs(x = "Half hour (start)",
y = "Mean halfhourly generation (GW)",
caption = ngEsoCap)
t <- ngeso_dt[, .(N = .N,
minTime = min(dv_hh_start_hms),
maxTime = max(dv_hh_start_hms),
meanGen_GW = mean(GENERATION),
meanCI = mean(CARBON_INTENSITY)), keyby = .(dv_peakPeriod)]
makeFlexTable(t, cap = "Half-hourly summaries by peak period (all data)")
dv_peakPeriod | N | minTime | maxTime | meanGen_GW | meanCI |
|---|---|---|---|---|---|
Early morning (00:00 - 06:00) | 65,100 | 0 | 19,800 | 27,543 | 293 |
Morning peak (06:00 - 09:00) | 32,550 | 21,600 | 30,600 | 35,200 | 323 |
Day time (09:00 - 16:00) | 75,950 | 32,400 | 55,800 | 39,088 | 322 |
Evening peak (16:00 - 20:00) | 43,395 | 57,600 | 70,200 | 39,790 | 337 |
Late evening (20:00 - 00:00) | 43,392 | 72,000 | 84,600 | 32,477 | 317 |
Check the response code. This seems to generate errors (sometimes).
# test
url <- "https://api.octopus.energy/v1/products"
message("Getting: ", url)
## Getting: https://api.octopus.energy/v1/products
resp <- httr::GET(url)
message("Status code: ", resp$status_code)
## Status code: 200
df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
makeFlexTable(head(df$results), cap = "Example products list (first 6 rows)")
code | direction | full_name | display_name | description | is_variable | is_green | is_tracker | is_prepay | is_business | is_restricted | term | available_from | available_to | links | brand |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
AGILE-FLEX-22-11-25 | IMPORT | Agile Octopus November 2022 v1 | Agile Octopus | With Agile Octopus, you get access to half-hourly energy prices, tied to wholesale prices and updated daily. The unit rate is capped at 100p/kWh (including VAT). | TRUE | TRUE | FALSE | FALSE | FALSE | FALSE | 2022-11-25T00:00:00Z | [[data.frame]] | OCTOPUS_ENERGY | ||
AGILE-FLEX-BB-23-02-08 | IMPORT | Agile Octopus February 2023 v1 | Agile Octopus | With Agile Octopus, you get access to half-hourly energy prices, tied to wholesale prices and updated daily. The unit rate is capped at 100p/kWh (including VAT). | TRUE | TRUE | FALSE | FALSE | FALSE | FALSE | 2023-02-08T00:00:00Z | [[data.frame]] | BULB | ||
AGILE-OUTGOING-19-05-13 | EXPORT | Agile Outgoing Octopus May 2019 | Agile Outgoing Octopus | Outgoing Octopus Agile rate pays you for all your exported energy based on the day-ahead wholesale rate. | TRUE | TRUE | FALSE | FALSE | FALSE | FALSE | 12 | 2018-01-01T00:00:00Z | [[data.frame]] | OCTOPUS_ENERGY | |
AGILE-OUTGOING-BB-23-02-28 | EXPORT | Agile Outgoing Octopus February 2023 v1 | Agile Outgoing Octopus | Outgoing Octopus Agile rate pays you for all your exported energy based on the day-ahead wholesale rate. | TRUE | TRUE | FALSE | FALSE | FALSE | FALSE | 12 | 2023-02-27T00:00:00Z | [[data.frame]] | BULB | |
COOP-LOYAL-FIX-12M-23-10-20 | IMPORT | Co-op Loyal 12M Fixed October 2023 v1 | Co-op Loyal 12M Fixed | This tariff features 100% renewable electricity and fixes your unit rates and standing charge for 12 months. This tariff has £75/fuel exit fees. | FALSE | FALSE | FALSE | FALSE | FALSE | FALSE | 12 | 2023-10-20T00:00:00+01:00 | [[data.frame]] | COOP_ENERGY | |
COOP-PP-VAR-20-04-01 | IMPORT | Co-op Key and Card | Co-op Key and Card | Non-smart prepayment tariff | TRUE | FALSE | FALSE | TRUE | FALSE | FALSE | 2020-03-03T00:00:00Z | [[data.frame]] | COOP_ENERGY |
url <- paste0("https://api.octopus.energy/v1/accounts/", apiParams$accountNo , "/")
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
props <- data.table::as.data.table(df$properties)
makeFlexTable(head(props[, .(town, county,
electricity_meter_points,gas_meter_points)]), cap = "Properties linked to this account (non-disclosive data)")
town | county | electricity_meter_points | gas_meter_points |
|---|---|---|---|
BRADFORD-ON-AVON | WILTSHIRE | [[data.frame]] | [[data.frame]] |
FRAMLINGHAM | [[data.frame]] | [[data.frame]] |
List the electricity meters by MPAN. There should be two - import & export…
# this is a list of n mpans
length(props$electricity_meter_points)
## [1] 2
message("n MPANS listed: ", length(df$properties$electricity_meter_points))
## n MPANS listed: 2
for(n in 1:length(df$properties$electricity_meter_points)){
print(props$electricity_meter_points[n])
}
## [[1]]
## mpan profile_class consumption_standard
## 1 2000006198482 1 2549
## meters
## 1 L78C64517, 01, STANDARD, TRUE
## agreements
## 1 E-1R-SUPER-GREEN-12M-20-09-22-H, E-1R-SUPER-GREEN-12M-20-09-22-H, 2020-11-01T00:00:00Z, 2020-11-21T00:00:00Z, 2020-11-01T00:00:00Z, 2021-06-30T00:00:00+01:00
## is_export
## 1 FALSE
##
## [[1]]
## mpan profile_class consumption_standard
## 1 1050002522164 8 1953
## 2 1050001805886 1 3374
## meters
## 1 19L3027004, 1, STANDARD, TRUE
## 2 19L3027004, 1, STANDARD, TRUE
## agreements
## 1 E-1R-AGILE-OUTGOING-19-05-13-A, E-1R-AGILE-OUTGOING-19-05-13-A, 2022-12-14T00:00:00Z, 2023-01-12T00:00:00Z, 2023-01-12T00:00:00Z, 2024-01-12T00:00:00Z
## 2 E-1R-SUPER-GREEN-12M-20-09-22-A, E-1R-LOYAL-FIX-12M-21-10-07-A, E-1R-VAR-22-10-01-A, E-1R-VAR-22-10-01-A, E-1R-VAR-22-11-01-A, 2021-07-01T00:00:00+01:00, 2021-11-21T00:00:00Z, 2022-11-21T00:00:00Z, 2022-12-21T00:00:00Z, 2023-04-01T00:00:00+01:00, 2021-11-21T00:00:00Z, 2022-11-21T00:00:00Z, 2023-04-01T00:00:00+01:00, 2022-12-21T00:00:00Z, NA
## is_export
## 1 TRUE
## 2 FALSE
List the gas meters by MPRN. There should be only one…
length(df$properties$gas_meter_points)
## [1] 2
message("n MPRNS listed: ", length(df$properties$gas_meter_points))
## n MPRNS listed: 2
for(n in 1:length(df$properties$gas_meter_points)){
print(df$properties$gas_meter_points[n])
}
## [[1]]
## mprn consumption_standard meters
## 1 4256845702 8586 G4A01559730801
## agreements
## 1 G-1R-SUPER-GREEN-12M-20-09-22-H, G-1R-SUPER-GREEN-12M-20-09-22-H, 2020-11-01T00:00:00Z, 2020-11-21T00:00:00Z, 2020-11-01T00:00:00Z, 2021-06-30T00:00:00+01:00
##
## [[1]]
## mprn consumption_standard meters
## 1 7825700304 15649 E6S12725512161, E6S17944211961, NOTINSTALLED
## agreements
## 1 G-1R-LOYAL-FIX-12M-21-10-07-A, G-1R-VAR-22-04-02-A, G-1R-VAR-22-10-01-A, G-1R-VAR-22-11-01-A, 2021-10-04T00:00:00+01:00, 2022-10-04T00:00:00+01:00, 2022-12-21T00:00:00Z, 2023-04-01T00:00:00+01:00, 2022-10-04T00:00:00+01:00, 2023-04-01T00:00:00+01:00, 2022-12-21T00:00:00Z, NA
use)See: https://www.guylipman.com/octopus/api_guide.html#s3
Start data extraction from 1st Jan 2022 as we had smart meter (re)installed in January.
Check missing dates and adjust “&page_size=100000” if required
url <- paste0("https://api.octopus.energy/v1/electricity-meter-points/",
apiParams$elec_import_mpan , "/",
"meters/",
apiParams$elec_import_serial, "/",
"consumption/",
"?period_from=2022-01-01T00:00Z",
"&page_size=100000") # make sure this is large enough!
# get data via httr ----
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE) # creates a df of which 'results' = the data
## No encoding supplied: defaulting to UTF-8.
elecCons_dt <- data.table::as.data.table(df$results) # convert to dt
# derived variables ----
elecCons_dt <- addDerivedVariables(elecCons_dt, source = "octopus")
maxTime <- max(elecCons_dt$dv_start)
hoursAgo <- now() - maxTime
# meter is SMETS2
elecCons_dt[, consumption_kWh := consumption] # for clarity - see https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter
message("# Check start and end dates")
## # Check start and end dates
summary(elecCons_dt$dv_start)
## Min. 1st Qu. Median
## "2022-06-05 04:00:00" "2022-10-13 08:52:30" "2023-02-20 13:45:00"
## Mean 3rd Qu. Max.
## "2023-02-20 13:45:00" "2023-06-30 18:37:30" "2023-11-07 23:30:00"
The data used here is up to 2023-11-07 23:30:00, which is 19.7 hours ago. In general the Octopus API seems to have data up to midnight last night.
Figure 5.1 shows half-hourly electricity import (‘consumption’) for the current year. Spot the power cuts…
To do: mark weekends somehow
ggplot2::ggplot(elecCons_dt, aes(x = dv_hh_start_date, y = dv_hh_start_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Electricity import (kWh)") +
labs(x = "Date",
y = "Half-hour")
Figure 5.1: Half hourly electricity import (current year)
Repeat but with just the last 14 days of data - useful for checking recent appliance use and offspring effects. We do quite a lot of batch cooking on Sunday nights…
Check this really is the last 14 days - there may be data errors
today <- lubridate::today()
plotDT <- elecCons_dt[dv_hh_start_date >= max(dv_hh_start_date) - 14]
p <- ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = dv_hh_start_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Electricity import (kWh)") +
scale_x_date(date_breaks = "1 day", date_labels = "%a %b %d") +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Date",
y = "Half-hour")
plotly::ggplotly(p)
Figure 5.2: Half hourly electricity import (current year, last 14 days)
plotDT[, dow := lubridate::wday(dv_hh_start_date, label = TRUE)]
recent_dt <- plotDT[, .(sum_kWh_elec = sum(consumption_kWh)), keyby = .(dv_hh_start_date, dow, dv_peakPeriod)]
daily_totals <- plotDT[, .(Total = sum(consumption_kWh)), keyby = .(dv_hh_start_date, dow)]
t <- dcast(recent_dt, dv_hh_start_date + dow ~ dv_peakPeriod, val.var = sum_kWh_elec)
## Using 'sum_kWh_elec' as value column. Use 'value.var' to override
# add totals
t <- t[daily_totals]
makeFlexTable(t, digits = 2,
cap = "Recent electricity use")
dv_hh_start_date | dow | Early morning (00:00 - 06:00) | Morning peak (06:00 - 09:00) | Day time (09:00 - 16:00) | Evening peak (16:00 - 20:00) | Late evening (20:00 - 00:00) | Total |
|---|---|---|---|---|---|---|---|
2023-10-24 | Tue | 1.36 | 1.12 | 2.10 | 7.20 | 2.12 | 13.89 |
2023-10-25 | Wed | 1.28 | 1.43 | 3.10 | 2.23 | 1.34 | 9.37 |
2023-10-26 | Thu | 1.30 | 1.40 | 1.70 | 2.33 | 2.25 | 8.98 |
2023-10-27 | Fri | 1.56 | 1.44 | 1.48 | 2.26 | 1.51 | 8.25 |
2023-10-28 | Sat | 1.22 | 0.73 | 1.37 | 3.62 | 2.17 | 9.11 |
2023-10-29 | Sun | 1.26 | 0.66 | 3.37 | 5.24 | 2.46 | 12.99 |
2023-10-30 | Mon | 1.48 | 1.20 | 3.54 | 3.02 | 3.86 | 13.09 |
2023-10-31 | Tue | 1.45 | 1.01 | 1.61 | 2.59 | 1.65 | 8.31 |
2023-11-01 | Wed | 1.20 | 0.84 | 3.77 | 3.61 | 2.48 | 11.91 |
2023-11-02 | Thu | 1.14 | 1.06 | 2.67 | 3.16 | 1.32 | 9.35 |
2023-11-03 | Fri | 1.51 | 0.92 | 1.03 | 3.15 | 2.91 | 9.52 |
2023-11-04 | Sat | 1.25 | 0.86 | 2.40 | 5.65 | 2.22 | 12.37 |
2023-11-05 | Sun | 1.65 | 0.81 | 0.89 | 2.69 | 1.87 | 7.91 |
2023-11-06 | Mon | 1.11 | 1.03 | 2.64 | 3.21 | 2.48 | 10.46 |
2023-11-07 | Tue | 1.12 | 1.12 | 3.63 | 2.24 | 1.24 | 9.35 |
daily_totals[, fuel := "elec"] # for future use
p <- ggplot2::ggplot(recent_dt, aes(x = dv_hh_start_date,
y = sum_kWh_elec, fill = dv_peakPeriod)) +
geom_col(position = "stack") +
scale_fill_viridis_d(name = "Time of day") +
scale_x_date(date_breaks = "1 day", date_labels = "%a %b %d") +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Date",
y = "Electricity kWh")
plotly::ggplotly(p)
Figure 5.2: Half hourly electricity import (current year, last 14 days)
plotDT <- elecCons_dt[dv_hh_start_date > Sys.Date()-14, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_hh_start_date)]
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = sum_kWh)) +
geom_col() +
#facet_grid(dv_peakPeriod ~ .) +
labs(x = "Date",
y = "Daily kWh")
Figure 5.3: Mean half-hourly electricity import per day (kWh, last 2 weeks)
Figure 5.4 shows the total daily kWh import with a smoothed curve for each day as shown.
elecCons_dt[, dv_month := lubridate::month(dv_hh_start_date, label = TRUE)]
elecCons_dt[, dv_year := lubridate::year(dv_hh_start_date)]
elecCons_dt[, dv_yday := lubridate::yday(dv_hh_start_date)]
plotDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_yday, dv_month, dv_year)]
makeDailyPlotByYear <- function(dt, yVar){
# expects a plotDT
p <- ggplot2::ggplot(plotDT, aes(x = dv_yday, y = get(yVar))) +
geom_smooth(aes(linetype = as.factor(dv_year))) +
geom_point(aes(shape = as.factor(dv_year),
colour = dv_month,
alpha = dv_year)) +
#facet_grid(dv_year ~ .) +
theme(legend.position = "bottom") +
scale_colour_viridis_d(name = "Month") +
scale_linetype_discrete(name = "Year") +
scale_alpha_continuous(name = "Year") +
scale_shape_discrete(name = "Year") +
scale_alpha_continuous(guide = "none") +
#scale_alpha_discrete(name = "Year") +
labs(x = "Day of the year",
y = "Daily kWh")
return(p)
}
makeDailyPlotByYear(plotDT, yVar = "sum_kWh")
## Scale for alpha is already present.
## Adding another scale for alpha, which will replace the existing scale.
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.4: Mean half-hourly electricity import per day (kWh, current year)
Figure 5.5 shows the mean daily kWh import with a smoothed curve for each period as defined below. The periods do not have the same number of half-hours so we use the mean as a comparator.
Early morning is effectively our baseload.
# check periods
t <- elecCons_dt[, .(min = min(dv_hh_start_hms),
max = max(dv_hh_start_hms)),
keyby = .(dv_peakPeriod)]
t
## dv_peakPeriod min max
## 1: Early morning (00:00 - 06:00) 00:00:00 05:30:00
## 2: Morning peak (06:00 - 09:00) 06:00:00 08:30:00
## 3: Day time (09:00 - 16:00) 09:00:00 15:30:00
## 4: Evening peak (16:00 - 20:00) 16:00:00 19:30:00
## 5: Late evening (20:00 - 00:00) 20:00:00 23:30:00
t <- elecCons_dt[, .(mean_kWh = mean(consumption_kWh),
sd_kWh = sd(consumption_kWh),
min_kWh = min(consumption_kWh),
max_kWh = max(consumption_kWh)),
keyby = .(dv_peakPeriod)]
makeFlexTable(t, digits = 3, cap = "Summary stats (all half-hourly data)")
dv_peakPeriod | mean_kWh | sd_kWh | min_kWh | max_kWh |
|---|---|---|---|---|
Early morning (00:00 - 06:00) | 0.124 | 0.064 | 0.044 | 1.038 |
Morning peak (06:00 - 09:00) | 0.177 | 0.131 | 0.000 | 1.273 |
Day time (09:00 - 16:00) | 0.129 | 0.180 | 0.000 | 1.994 |
Evening peak (16:00 - 20:00) | 0.327 | 0.285 | 0.000 | 1.720 |
Late evening (20:00 - 00:00) | 0.274 | 0.205 | 0.000 | 1.900 |
plotDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_hh_start_date, dv_peakPeriod)]
totalDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_hh_start_date)]
totalDT[, dv_peakPeriod := "All periods"]
plotDT <- rbind(plotDT, totalDT)
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = mean_kWh,
colour = dv_peakPeriod)) +
geom_line() +
geom_smooth() +
#facet_grid(dv_peakPeriod ~ .) +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Peak period") +
labs(x = "Date",
y = "Mean kWh per half-hour")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.5: Mean half-hourly electricity import (kWh) by peak period (current year)
What’s happening here?
Day timeevening peakearly morning baseload) and also swapped a built-in fridge-freezer for a stand-alone fridge (& freezer)morning peak). Well, kind of. Might affect the gas more.Not Octopus data - irregular reads by us. Load from .xslsx.
To do: update to use octopus Agile Outgoing data
Figure 5.6 compares our PV generation with overall grid import. These are mean values. We are occasionally exporting (Table ?? has zeros) but we do not (currently) have access to the export data.
library(openxlsx)
f <- "~/Dropbox/Home/houses/whiteHorseMews/2WhiteHorseMewsRunningCosts.xlsx"
pvGen <- openxlsx::read.xlsx(f,
sheet = "PV",
detectDates = TRUE)
pvGen_DT <- data.table::as.data.table(pvGen)
pvGen_DT[, meankWhGen := diff/n.days]
pvGen_DT[, dv_hh_start_date := lubridate::as_datetime(Date)]
pvGen_DT[, month := lubridate::month(dv_hh_start_date)]
monthlyPvGen_DT <- pvGen_DT[, .(hh_meankWh = mean(meankWhGen)/48 # input value is mean per day
),
keyby = .(month = lubridate::month(dv_hh_start_date),
year = lubridate::year(dv_hh_start_date))
]
monthlyPvGen_DT[, value := "PV generation"]
monthlyElec_DT <- elecCons_dt[, .(hh_meankWh = mean(consumption_kWh)),
keyby = .(month = lubridate::month(dv_hh_start_date),
year = lubridate::year(dv_hh_start_date))]
monthlyElec_DT[, value := "Grid import"]
plotDT <- rbind(monthlyPvGen_DT, monthlyElec_DT)
ggplot2::ggplot(plotDT, aes(x = lubridate::month(month, label = TRUE), y = hh_meankWh,
colour = value, group = value)) +
geom_line()+
scale_color_discrete(name = "Type")+
facet_grid(year ~ .) +
labs(x = "Month")
## Warning: Removed 1 row containing missing values (`geom_line()`).
Figure 5.6: Compare grid import and PV gen
To do: model value of PV gen if we were to use all of it.
This will be a new MPAN but specified as export - although the url will still say consumption. We do not have this even though the PV is exporting on (some) days.
It may be that we only get this data if we sign up for the export tariff.
In theory our emissions from electricity use are zero because we are on a renewable-only tariff. But life is not so simple. We don’t have a private wire to a wind turbine so the electrons we import (stick with it) are as averagely green as all the rest.
We also con’t avoid the ‘Well To Tank’ emissions and those associated with transmission losses.
To further complicate things there are at least two different ways to estimate our emissions.
Does it matter? you cry. Well it might. If we’ve been able to ‘flex’ our usage in line with @theBakingForecast then who knows, maybe we’ll be concentrating our usage in times when the grid is actually drawing on more renewables.
So let’s take a look. We’ll do both the BEIS-based and NG-ESO based calculations to see. For now we’ll ignore the WTT and the T&D losses to keep the results comparable. We’ll come back to that later.
rmdParams$BEIS_elec_ci <- 0.21233
For the BEIS method, we’ll have to use the 2021 emissions factor as the 2022 value is not yet available.
For 2021 this is: 0.21233 Kg CO2e/kWh
For the NG-ESO method we use the NG-ESO half-hourly carbon intensity data that match to the half-hours in our electricity use dataset.
Mean half-hourly carbon intensity from the NG-ESO data for the data period was NA Kg CO2e/kWh. If this is substantially different to the BEIS 2021 value of 0.2123 Kg CO2e/kWh, we would expect emissions estimates using the NG-ESO factor to differ.
# merge to usage data
setkey(ngeso_dt, dv_start)
setkey(elecCons_dt, dv_start)
elecCons_dt <- ngeso_dt[, .(dv_start, CARBON_INTENSITY, LOW_CARBON_perc, RENEWABLE_perc)][elecCons_dt] # keeps match to our electricity use
# there will be NAs if some datetimes are missing from ngeso_dt
For context, Figure 5.7 summarises the mean half-hourly carbon intensity by month for the data period. We can clearly see that February 2022 was a very low carbon month… in fact it was a very windy month with 3 named storms.
elecCons_dt[, dv_month := lubridate::month(dv_hh_start_date, label = TRUE)]
ggplot2::ggplot(elecCons_dt, aes(x = dv_month, y = CARBON_INTENSITY)) +
geom_violin(draw_quantiles = c(0.5)) +
facet_grid(dv_year ~ .) +
#geom_jitter() +
#geom_boxplot() +
labs(x = "Month",
y = "Half-hourly carbon intensity",
caption = "Median drawn")
Figure 5.7: Monthly mean carbon intensity for the data period by month (NG-ESO data)
Figure 5.8 shows our half-hourly electricity kWh use vs halfhourly carbon intensity. Ideally we want a negative correlation showing that we use the most electricity when it is ‘greenest’ (carbon intensity is lowest). Doesn’t look too good, aye?
ggplot2::ggplot(elecCons_dt, aes(x = CARBON_INTENSITY, y = consumption_kWh, colour = RENEWABLE_perc)) +
geom_point() +
facet_wrap(. ~ dv_peakPeriod) +
geom_smooth() +
scale_color_continuous(name = "% renewables", low = "red", high = "green") +
theme(legend.position = "bottom") +
labs(y = "Halfhourly electricity kWh")
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
Figure 5.8: Own half-hourly electricity kWh vs NG-ESO halfhourly carbon intensity
What if we visualise using a box plot according to carbon intensity decile? So this means we divide the carbon intensity values into 10 equal groups - deciles. This is Figure 5.9. Doesn’t look too good either - our median usage (the horizontal bar in the boxes) seems to trend slightly upwards as we move to higher carbon intensity deciles.
# this will generate NAs if the CI data is missing for some of the (most recent) dateTimes
elecCons_dt[, CI_deciles := cut_number(CARBON_INTENSITY, n = 10)]
ggplot2::ggplot(elecCons_dt, aes(x = CI_deciles, y = consumption_kWh)) +
geom_boxplot() +
labs(y = "Halfhourly electricity kWh",
x = "Carbon intensity decile")
Figure 5.9: Boxing clever
So what if we just add up all our electricity kWh by carbon intensity decile? Do we use more low carbon kWh? This is Figure 5.10. Nah. The bakingforecast isn’t going to like us…
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh, na.rm = TRUE),
meankWh = mean(consumption_kWh, na.rm = TRUE)),
keyby = .(CI_deciles)]
ggplot2::ggplot(t, aes(x = CI_deciles, y = sumkWh)) +
geom_col() +
labs(y = "Sum kWh",
x = "Carbon intensity decile")
Figure 5.10: Sum of electricity kWh by carbon intensity decile
Out of interest, do our emissions values look very different if we apply the BEIS 2021 annual factor to our total electricity kWh to date compared to applying the NG-ESO half-hourly values?
elecCons_dt[, KgCO2_ngeso := consumption_kWh * (CARBON_INTENSITY/1000)] # convert to kg
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh),
sumKgCO2_ngeso = sum(KgCO2_ngeso, na.rm = TRUE)), keyby = .(Year = as.factor(dv_year))]
t[, sumKgCO2_beis := sumkWh * rmdParams$BEIS_elec_ci]
makeFlexTable(t, cap = "Comparing emissions estimation methods using electricity kWh to date")
Year | sumkWh | sumKgCO2_ngeso | sumKgCO2_beis |
|---|---|---|---|
2022 | 2,144 | 405 | 455 |
2023 | 2,627 | 415 | 558 |
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh),
sumKgCO2_ngeso = sum(KgCO2_ngeso, na.rm = TRUE)),
keyby = .(dv_month, dv_year)]
t[, sumKgCO2_beis := sumkWh * rmdParams$BEIS_elec_ci]
plotDT <- melt(t, id.vars = c("dv_month", "dv_year"))
ggplot2::ggplot(plotDT[ variable != "sumkWh",], aes(x = dv_month, y = value, fill = variable)) +
geom_col(position = "dodge") +
facet_grid(dv_year ~ .) +
scale_color_discrete(name = "Method") +
labs(y = "Kg CO2",
x = "Month")
As we’d expect from the comparison of the values above, Table 5.4 suggests that it does. In fact our ‘in use’ NG-ESO based emissions are 67.29, 81.05, 78.66, 77.33, 75.3, 87.38, 81.97, 99.48, 67.58, 105.82, 78.58, 100.65, 76.1, 76.93, 68.34, 83.46, 58.23, 73.77 % of our BEIS-based emissions depending on the month in question.
If we compare the monthly values we can see the biggest difference was in February, a month we have already identified as being more ‘low carbon’ (see Figure 5.7).
Analyse costs using:
The latter are slightly different from the assumed to be at UK price cap: £0.34 / kWh & £0.46 ( see Ofgem)
prices <- readxl::read_xlsx(here::here("data", "prices.xlsx"))
pricesDT <- data.table::as.data.table(prices)
Yes, I know I can extract our exact tariff from the octopus API…
daily_elec <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh,
na.rm = TRUE), # beware missing (N/A) may decrease sum
nObs = .N), keyby = .(dv_hh_start_date, dv_yday, dv_month, dv_year)]
# extract from pricesDT
# must be an easier way
daily_elec[, kwh_p := ifelse(dv_hh_start_date < lubridate::as_date("2022-11-21"),
pricesDT[fuel == "elec_imp" & component == "kWh" &
dateEnd == lubridate::as_datetime("2022-11-21"), price],
pricesDT[fuel == "elec_imp" & component == "kWh" &
dateStart == lubridate::as_datetime("2022-11-22"), price])]
daily_elec[, sc_p := ifelse(dv_hh_start_date < lubridate::as_date("2022-11-21"),
pricesDT[fuel == "elec_imp" & component == "sc" &
dateEnd == lubridate::as_datetime("2022-11-21"), price],
pricesDT[fuel == "elec_imp" & component == "sc" &
dateStart == lubridate::as_datetime("2022-11-22"), price])]
daily_elec[, cost := ((sum_kWh * kwh_p) + sc_p)]
daily_elec[, month := lubridate::month(dv_hh_start_date, label = TRUE)]
ggplot2::ggplot(daily_elec, aes(x = dv_hh_start_date, y = cost)) +
geom_point(aes(colour = month)) +
geom_smooth() +
geom_vline(xintercept = lubridate::as_date("2022-11-21")) +
labs(y = "Electricity daily cost £",
caption = "Tariff change/price cap/EPG 2022 date shown\nSmoothed within month")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.11: Daily electricity costs
lastWeek <- max(daily_elec$dv_hh_start_date) - 7
makeFlexTable(daily_elec[dv_hh_start_date > lastWeek, .(dv_hh_start_date,
day = lubridate::wday(dv_hh_start_date, label = TRUE),
sum_kWh, nObs, cost)], digits = 2,
cap = "Recent daily gas cost")
dv_hh_start_date | day | sum_kWh | nObs | cost |
|---|---|---|---|---|
2023-11-01 | Wed | 11.91 | 48 | 4.25 |
2023-11-02 | Thu | 9.35 | 48 | 3.43 |
2023-11-03 | Fri | 9.52 | 48 | 3.48 |
2023-11-04 | Sat | 12.37 | 48 | 4.40 |
2023-11-05 | Sun | 7.91 | 48 | 2.96 |
2023-11-06 | Mon | 10.46 | 48 | 3.79 |
2023-11-07 | Tue | 9.34 | 48 | 3.42 |
#makeDailyPlotByYear(daily_elec, yVar = "cost") # should work but doesn't
daily_elec[, month_floor := lubridate::floor_date(dv_hh_start_date, "months")]
monthly_elec <- daily_elec[, .(sum_kWh = sum(sum_kWh),
cost = sum(cost)),
keyby = .(month_floor)]
ggplot2::ggplot(monthly_elec, aes(x = month_floor, y = cost)) +
geom_col() +
geom_vline(xintercept = lubridate::as_date("2022-11-01")) +
labs(y = "Monthly cost £",
x = "Month",
caption = "Tariff change/price cap/EPG 2022 date shown\nBeware incomplete months")
Figure 5.12: Monthly electricity costs
message("Projected annual elec total kWh")
## Projected annual elec total kWh
projAnnual_elec_kWh <- mean(daily_elec$sum_kWh)*365
projAnnual_elec_kWh
## [1] 3342.04
message("########")
## ########
message("# Prices to 21st Nov 2022")
## # Prices to 21st Nov 2022
# get_price <- function(dt, fuel, component, dateEnd){
# p <- dt[fuel == fuel &
# component == component &
# dateEnd == dateEnd,
# price]
# return(p)
# }
kWh_p <- pricesDT[fuel == "elec_imp" & component == "kWh" &
dateEnd == lubridate::as_datetime("2022-11-21"), price]
sc_p <- pricesDT[fuel == "elec_imp" & component == "sc" &
dateEnd == lubridate::as_datetime("2022-11-21"), price]
message("Projected annual elec cost @ £ ", kWh_p,
" per kWh & standing charge at £ ", sc_p,
" per day")
## Projected annual elec cost @ £ 0.2408 per kWh & standing charge at £ 0.2401 per day
projAannualCost <- (projAnnual_elec_kWh*kWh_p)+(365*sc_p) # standing charge
projAannualCost
## [1] 892.3998
message("Mean projected monthly cost: £")
## Mean projected monthly cost: £
projAannualCost/12
## [1] 74.36665
message("########")
## ########
message("# From to 22nd Nov 2022 - price cap")
## # From to 22nd Nov 2022 - price cap
kWh_p <- pricesDT[fuel == "elec_imp" & component == "kWh" &
dateStart == lubridate::as_datetime("2022-11-22"), price]
sc_p <- pricesDT[fuel == "elec_imp" & component == "sc" &
dateStart == lubridate::as_datetime("2022-11-22"), price]
message("Projected annual elec cost @ £ ", kWh_p,
" per kWh & standing charge at £ ", sc_p,
" per day")
## Projected annual elec cost @ £ 0.3236 per kWh & standing charge at £ 0.4001 per day
annualCost_elecCapped <- (projAnnual_elec_kWh*kWh_p)+(365*sc_p) # standing charge
annualCost_elecCapped
## [1] 1227.521
message("Projected mean monthly £ under price cap")
## Projected mean monthly £ under price cap
monthlyCost_capped <- annualCost_elecCapped/12
monthlyCost_capped
## [1] 102.2934
message("That's an increase of ", round(100*((annualCost_elecCapped-projAannualCost)/projAannualCost),2), " % points")
## That's an increase of 37.55 % points
NB - these willbe out of date from July 2023 as prices are changing
TO DO: fix £ analysis to use tariff from API
URL will be something like
Can we flex any of the above for ~ £3/kWh?
Instigated by UK National gird, implemented (in our case) by Octopus- https://twitter.com/SavingSessions
pd <- function(t){ # make a pretty date
strftime(t, "%a %d %b %Y")
}
pdt <- function(t){ # make a pretty date
strftime(t, "%a %d %b %Y %H:%M")
}
# re-usable plot function
make_kwhComparisonPlot <- function(dt, # half hourly smart meter data
startDateTime, # start
endDateTime, # end
timeLag = 10){ # n days to compare against
res <- list() # results holder
sessionDate <- lubridate::date(startDateTime) # you'll see why
if(lubridate::wday(startDateTime) < 6){
sessionDay <- "weekday" # you'll see why
} else {
sessionDay <- "weekend"
}
# 1 = Sunday, 7 = Saturday etc
dt[, ba_wday := lubridate::wday(dv_start, label = TRUE)]
dt[, ba_wd := ifelse(lubridate::wday(dv_start) > 1 & lubridate::wday(dv_start) < 7,
"weekday",
"weekend")]
# extract and aggregate baseline
# first get similar days (weekdays or weekends as needed)
similarDays <- dt[dv_hh_start_date != sessionDate & ba_wd == sessionDay &
dv_hh_start_date < sessionDate]
uniqueN(similarDays$dv_hh_start_date)
dates <- tail(similarDays[,(nObs = .N), keyby = .(dv_hh_start_date)], timeLag) # get the most recent n = timeLag days
datesToGet <- dates[, dv_hh_start_date]
daysWeWant <- similarDays[dv_hh_start_date %in% datesToGet]
uniqueN(daysWeWant$dv_hh_start_date)
baseline_meanHalfHourlyDT <- daysWeWant[,
.(mean_kWh_elec = mean(consumption_kWh),
minDate = min(dv_hh_start_date),
maxDate = max(dv_hh_start_date)),
keyby = .(hms = dv_hh_start_hms)]
res$nCompDays <- uniqueN(daysWeWant$dv_hh_start_date)
baseline_meanHalfHourlyDT[, legend_lab := paste0("Comparison mean")]
# extract session day
session_HalfHourlyDT <- dt[dv_hh_start_date == lubridate::date(startDateTime), # the session day
.(mean_kWh_elec = mean(consumption_kWh),
minDate = min(dv_hh_start_date),
maxDate = max(dv_hh_start_date)),
keyby = .(hms = dv_hh_start_hms)]
session_HalfHourlyDT[, legend_lab := "Saving session"]
plotDT <- rbind(baseline_meanHalfHourlyDT, session_HalfHourlyDT)
plotDT[, adjusted_hms := hms::as_hms(hms + (15*60))] # plots points in centre of half-hour period for clarity
periodAlpha <- 0.3 # shaded rects on plots
periodFill <- "grey50"
ymax <- max(plotDT$mean_kWh_elec)
ymin <- min(plotDT$mean_kWh_elec)
xmin <- hms::as_hms(lubridate::as_datetime(startDateTime))
xmax <- hms::as_hms(lubridate::as_datetime(endDateTime)) + 30*60 # to allow for the start time
# make the table before we muck about for the plot
t <- plotDT[hms >= xmin & hms < xmax]
wt <- dcast(t[, .(hms, wday, mean_kWh_elec, legend_lab)], hms ~ legend_lab, value.var = "mean_kWh_elec")
wt[, kwh_diff := `Saving session` - `Comparison mean`]
res$wt <- wt[, pc_diff := 100*(kwh_diff/`Comparison mean`)] # add to results holder for return
label <- paste0("Saving session (",
as.Date(startDateTime),
")"
)
plotDT[, legend_lab := ifelse(legend_lab == "Saving session",
label,
legend_lab)]
res$p <- ggplot2::ggplot(plotDT, aes(x = adjusted_hms, y = mean_kWh_elec,
colour = legend_lab)) +
geom_line() +
geom_point() +
annotate("rect", xmin = xmin,
xmax = xmax,
ymin = ymin, ymax = ymax,
alpha = periodAlpha, fill = periodFill) +
scale_color_manual(name = "Legend", values=c('grey', 'red')) + # this should always make comparison grey
theme(legend.position = "bottom") +
labs(x = "Time of day",
y = "Mean kWh per half-hour",
caption = paste0("Data: @OctopusEnergy\nPlot: @dataknut\nPoints centered in half-hours for clarity\nComparison = ", res$nCompDays," previous days (",
pd(min(daysWeWant$dv_hh_start_date)), " - ",
pd(max(daysWeWant$dv_hh_start_date)), ")"
)
)
return(res)
}
startDateTime <- "2022-11-15 17:00:00" # the half-hour it starts
endDateTime <- "2022-11-15 17:30:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.50
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Tue 15 Nov 2022 17:00 - 18:00
Rate: £2.5/kWh
Figure 5.13 shows how we did in the first #SavingSessions compared to the last 10 similar days.
NB: this may not be quite the same as the Octopus algorithm
res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Figure 5.13: First #SavingSessions - how did we do?
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
sumSavings <- totalSavingGBP
Total saving: 0.2066 kWh
Total prize: £ 0.52
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days)"),
digits = 3) # the table
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:00:00 | 0.359 | 0.386 | 0.027 | 7.401 |
17:30:00 | 0.226 | 0.406 | 0.180 | 79.646 |
What do we conclude?
For a complicated set of pre-Xmas prep reasons (and poor institutional memory) our electricity usage in the late afternoon was much higher than usual on the session day. As a result our actions only brought usage down to ‘average’ for this time on a Tuesday (so we won’t win any points). But as a result our usage was probably way lower than it otherwise would have been judging by the spikes before and after. So if it really had been a critical peak event, we would have been helping the system but not being rewarded…
startDateTime <- "2022-11-22 17:30:00" # the half-hour it starts
endDateTime <- "2022-11-22 18:00:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.25
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Tue 22 Nov 2022 17:30 - 18:30
Rate: £2.25/kWh
We missed this one but Figure 5.14 shows how we would have done in this #SavingSession compared to the last n similar days of the week - had we signed up.
if(as.Date(startDateTime) < Sys.Date()){ # octopus data only available for day before
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
data <- FALSE
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.14: Second #SavingSessions - how would we have done?
What do we conclude?
Total saving: -0.083 kWh
Total prize: £ -0.19
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:30:00 | 0.240 | 0.324 | 0.084 | 35.056 |
18:00:00 | 0.541 | 0.374 | -0.167 | -30.882 |
startDateTime <- "2022-11-30 17:30:00" # the half-hour it starts
endDateTime <- "2022-11-30 18:00:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.50
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Wed 30 Nov 2022 17:30 - 18:30
Rate: £2.5/kWh
Figure 5.15 shows how we did in the third #SavingSessions compared to the last n similar days of the week.
res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Figure 5.15: Third #SavingSessions - how did we do?
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
sumSavings <- totalSavingGBP + sumSavings
What do we conclude?
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:30:00 | 0.294 | 0.148 | -0.146 | -49.694 |
18:00:00 | 0.474 | 0.148 | -0.326 | -68.783 |
startDateTime <- "2022-12-01 17:00:00" # the half-hour it starts
endDateTime <- "2022-12-01 17:30:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.50
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Thu 01 Dec 2022 17:00 - 18:00
Rate: £2.5/kWh
Figure 5.16 shows how we did in this #SavingSession compared to the last n similar days of the week.
res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Figure 5.16: Fourth #SavingSessions - how did we do?
What do we conclude?
Total saving: -0.2558 kWh
Total prize: £ -0.64
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:00:00 | 0.355 | 0.179 | -0.176 | -49.563 |
17:30:00 | 0.290 | 0.210 | -0.080 | -27.561 |
startDateTime <- "2022-12-12 17:00:00" # the half-hour it starts
endDateTime <- "2022-12-12 17:30:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.50
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Mon 12 Dec 2022 17:00 - 18:00
Rate: £2.5/kWh
Figure 5.17 shows how we did in this #SavingSession compared to the last n similar days of the week.
res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Figure 5.17: Fifth #SavingSessions - how did we do?
What do we conclude?
Total saving: -0.4243 kWh
Total prize: £ -1.06
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:00:00 | 0.357 | 0.114 | -0.243 | -68.085 |
17:30:00 | 0.279 | 0.098 | -0.181 | -64.887 |
startDateTime <- "2023-01-19 09:00:00" # the half-hour it starts
endDateTime <- "2023-01-19 09:30:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.50
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Thu 19 Jan 2023 09:00 - 10:00
Rate: £2.5/kWh
Figure 5.18 shows how we would have done had we entered for this #SavingSession compared to the last n similar days of the week.
res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Figure 5.18: Sixth #SavingSessions - how would we have done?
What do we conclude?
Total saving: -0.0748 kWh
Total prize: £ -0.19
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
09:00:00 | 0.225 | 0.161 | -0.064 | -28.540 |
09:30:00 | 0.204 | 0.194 | -0.010 | -5.134 |
startDateTime <- "2023-01-23 17:00:00" # the half-hour it starts
endDateTime <- "2023-01-23 17:30:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 3.37
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Mon 23 Jan 2023 17:00 - 18:00
Rate: £3.37/kWh
Figure 5.19 shows how we did in this #SavingSession compared to the last n similar days of the week.
if(as.Date(startDateTime) < Sys.Date()-1){ # data only available 24 hours later
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.19: Seventh #SavingSessions - how did we do?
What do we conclude?
Total saving: -0.2353 kWh
Total prize: £ -0.79
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:00:00 | 0.259 | 0.121 | -0.138 | -53.264 |
17:30:00 | 0.257 | 0.160 | -0.097 | -37.840 |
startDateTime <- "2023-01-24 16:30:00" # the half-hour it starts
endDateTime <- "2023-01-24 17:30:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 4
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Tue 24 Jan 2023 16:30 - 18:00
Rate: £4/kWh
Figure 5.20 shows how we did in this #SavingSession compared to the last n similar days of the week.
if(as.Date(startDateTime) < Sys.Date()){ # octopus data only available for day before
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
data <- FALSE
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.20: Eighth #SavingSessions - how did we do?
What do we conclude?
Total saving: -0.3641 kWh
Total prize: £ -1.46
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
16:30:00 | 0.272 | 0.135 | -0.137 | -50.386 |
17:00:00 | 0.250 | 0.151 | -0.099 | -39.576 |
17:30:00 | 0.248 | 0.120 | -0.128 | -51.632 |
startDateTime <- "2023-01-30 09:00:00" # the half-hour it starts
endDateTime <- "2023-01-30 09:30:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.25
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Mon 30 Jan 2023 09:00 - 10:00
Rate: £2.25/kWh
We missed this one but Figure 5.21 shows how we did in this #SavingSession compared to the last n similar days of the week - had we signed up.
if(as.Date(startDateTime) < Sys.Date()){ # octopus data only available for day before
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
data <- FALSE
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.21: Ninth #SavingSessions - how would we have done?
What do we conclude?
Total saving: 0.2505 kWh
Total prize: £ 0.56
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
09:00:00 | 0.151 | 0.130 | -0.021 | -13.964 |
09:30:00 | 0.182 | 0.454 | 0.272 | 148.904 |
startDateTime <- "2023-02-13 17:30:00" # the half-hour it starts
endDateTime <- "2023-02-13 18:00:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.25
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Mon 13 Feb 2023 17:30 - 18:30
Rate: £2.25/kWh
Figure 5.22 shows how we did in this #SavingSession compared to the last n similar days of the week.
if(as.Date(startDateTime) < Sys.Date()){ # octopus data only available for day before
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
data <- FALSE
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.22: Tenth #SavingSessions- how did we do?
What do we conclude?
Total saving: -0.3675 kWh
Total prize: £ -0.83
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:30:00 | 0.295 | 0.175 | -0.120 | -40.678 |
18:00:00 | 0.441 | 0.194 | -0.247 | -56.059 |
startDateTime <- "2023-02-21 17:30:00" # the half-hour it starts
endDateTime <- "2023-02-21 18:00:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.25
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Tue 21 Feb 2023 17:30 - 18:30
Rate: £2.25/kWh
Figure 5.23 shows how we would have done in this #SavingSession compared to the last n similar days of the week.
if(as.Date(startDateTime) < Sys.Date()){ # octopus data only available for day before
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
data <- FALSE
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.23: Eleventh #SavingSessions- how would we have done?
What do we conclude?
Total saving: -0.0594 kWh
Total prize: £ -0.13
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
17:30:00 | 0.236 | 0.210 | -0.026 | -10.866 |
18:00:00 | 0.329 | 0.295 | -0.034 | -10.280 |
Would have been better yesterday - power cut 14:00 - 21:00!
startDateTime <- "2023-03-15 18:30:00" # the half-hour it starts
endDateTime <- "2023-03-15 19:00:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 2.25
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Wed 15 Mar 2023 18:30 - 19:30
Rate: £2.25/kWh
Figure 5.24 shows how did in this #SavingSession compared to the last n similar days of the week.
if(as.Date(startDateTime) < Sys.Date()){ # octopus data only available for day before
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
data <- FALSE
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.24: Twelth #SavingSessions- how did we do?
What do we conclude?
Total saving: -0.4655 kWh
Total prize: £ -1.05
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
18:30:00 | 0.366 | 0.166 | -0.200 | -54.595 |
19:00:00 | 0.358 | 0.092 | -0.266 | -74.294 |
startDateTime <- "2023-03-23 18:30:00" # the half-hour it starts
endDateTime <- "2023-03-23 19:00:00" # the half-hour it ends (makes data selection easier)
gbp_rate <- 3.00
periodStr <- paste0(pdt(startDateTime), " - ", strftime(lubridate::as_datetime(endDateTime)+30*60, "%H:%M"))
From Thu 23 Mar 2023 18:30 - 19:30
Rate: £3/kWh
Figure 5.25 shows how did in this #SavingSession compared to the last n similar days of the week.
if(as.Date(startDateTime) < Sys.Date()){ # octopus data only available for day before
data <- TRUE
res <- make_kwhComparisonPlot(elecCons_dt,
startDateTime, endDateTime)
print(res$p) # the plot
totalSavingKWh <- sum(res$wt$kwh_diff)
totalSavingGBP <- round(sum(res$wt$kwh_diff) * gbp_rate,2)
} else {
data <- FALSE
theMsg <- "Data not yet available, try later :-("
message(theMsg)
totalSavingKWh <- theMsg
totalSavingGBP <- theMsg
}
Figure 5.25: Thirteenth #SavingSessions- how did we do?
What do we conclude?
Total saving: -0.3772 kWh
Total prize: £ -1.13
if(data){ # set above
makeFlexTable(res$wt[, hms:= as.factor(hms)],
cap = paste0("kWh comparisons (hms = half hour period start, comparison = ", res$nCompDays," previous days"),
digits = 3) # the table
} else {
message("Data not yet available, try later :-(")
}
hms | Comparison mean | Saving session | kwh_diff | pc_diff |
|---|---|---|---|---|
18:30:00 | 0.247 | 0.123 | -0.124 | -50.243 |
19:00:00 | 0.371 | 0.118 | -0.253 | -68.194 |
use)We need to convert the gas consumption from m3 to kWh - see https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter
gasM3TokWh <- 11.36
We use a multiplier of 11.36 kWh/m3 (https://www.theenergyshop.com/guides/how-to-convert-gas-units-to-kwh)
Check for missing dates and adjust “&page_size=100000” if needed
url <- paste0("https://api.octopus.energy/v1/gas-meter-points/",
apiParams$gas_mpan , "/",
"meters/",
apiParams$gas_serial, "/",
"consumption",
"?period_from=2022-01-01T00:00Z",
"&page_size=100000") # make sure is large enough
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
gasCons_dt <- data.table::as.data.table(df$results)
gasCons_dt <- addDerivedVariables(gasCons_dt, source = "octopus")
# gas 'consumption' is m3 - https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter
# convert to kWh
gasCons_dt[, consumption_m3 := consumption]
gasCons_dt[, consumption_kWh := consumption * gasM3TokWh]
message("# Check start and end dates")
## # Check start and end dates
summary(gasCons_dt$dv_start)
## Min. 1st Qu. Median
## "2022-05-07 06:30:00.000" "2022-09-14 11:22:30.000" "2023-01-30 17:15:00.000"
## Mean 3rd Qu. Max.
## "2023-02-06 23:39:56.375" "2023-06-29 05:37:30.000" "2023-11-07 23:00:00.000"
maxTime <- max(gasCons_dt$dv_start)
hoursAgo <- lubridate::now() - maxTime
Note that this data starts later as we finally got the original un-registered smart meter replaced in February 2022.
The data used here is up to 2023-11-07 23:00:00. In general the Octopus API seems to have data up to midnight last night.
Figure 5.26 shows half-hourly gas import (‘consumption’) for the current year. The power cuts are even easier to see here.
Gas boiler services:
ggplot2::ggplot(gasCons_dt, aes(x = dv_hh_start_date, y = dv_hh_start_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Gas import (kWh)") +
labs(x = "Date",
y = "Half-hour")
Figure 5.26: Half-hourly gas consumption (current year)
Repeat but with just the last 14 days of data - useful for checking recent appliance use and offspring effects.
Check this really is the last 2 weeks of gas use - there may be gas data errors
today <- lubridate::today()
plotDT <- gasCons_dt[dv_hh_start_date >= max(dv_hh_start_date) - 14] # last 14 days of data - might not be last 14 days in reality due to data errors
p <- ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = dv_hh_start_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Gas import (kWh)") +
scale_x_date(date_breaks = "1 day", date_labels = "%a %b %d") +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Date",
y = "Half-hour")
plotly::ggplotly(p)
Figure 5.27: Half hourly gas import (current year, last 14 days)
plotDT[, dow := lubridate::wday(dv_hh_start_date, label = TRUE)]
recent_dt <- plotDT[, .(sum_kWh_gas = sum(consumption_kWh)), keyby = .(dv_hh_start_date, dow, dv_peakPeriod)]
daily_totals <- plotDT[, .(sum_kWh_gas = sum(consumption_kWh)), keyby = .(dv_hh_start_date, dow)]
t <- dcast(recent_dt, dv_hh_start_date + dow ~ dv_peakPeriod, val.var = sum_kWh_gas)
## Using 'sum_kWh_gas' as value column. Use 'value.var' to override
# add totals
t <- t[daily_totals]
makeFlexTable(t, digits = 2,
cap = "Recent gas use")
dv_hh_start_date | dow | Early morning (00:00 - 06:00) | Morning peak (06:00 - 09:00) | Day time (09:00 - 16:00) | Evening peak (16:00 - 20:00) | Late evening (20:00 - 00:00) | sum_kWh_gas |
|---|---|---|---|---|---|---|---|
2023-10-24 | Tue | 15.87 | 5.00 | 8.87 | 12.77 | 0.00 | 42.51 |
2023-10-25 | Wed | 14.51 | 5.59 | 8.10 | 4.98 | 0.00 | 33.17 |
2023-10-26 | Thu | 16.32 | 10.79 | 8.82 | 8.09 | 0.86 | 44.88 |
2023-10-27 | Fri | 22.16 | 10.87 | 8.52 | 17.41 | 0.00 | 58.97 |
2023-10-28 | Sat | 14.67 | 10.92 | 8.52 | 7.35 | 0.00 | 41.45 |
2023-10-29 | Sun | 8.89 | 6.42 | 0.00 | 19.40 | 0.00 | 34.72 |
2023-10-30 | Mon | 8.88 | 19.40 | 0.00 | 20.89 | 0.20 | 49.38 |
2023-10-31 | Tue | 9.30 | 22.17 | 0.00 | 23.36 | 0.00 | 54.83 |
2023-11-01 | Wed | 9.53 | 14.25 | 0.00 | 15.70 | 0.00 | 39.48 |
2023-11-02 | Thu | 9.85 | 21.97 | 0.00 | 35.38 | 0.83 | 68.02 |
2023-11-03 | Fri | 7.94 | 19.94 | 0.00 | 27.97 | 0.00 | 55.85 |
2023-11-04 | Sat | 8.40 | 19.07 | 0.00 | 24.08 | 0.00 | 51.55 |
2023-11-05 | Sun | 8.72 | 14.78 | 0.00 | 27.04 | 0.22 | 50.76 |
2023-11-06 | Mon | 9.39 | 16.84 | 0.00 | 22.89 | 0.00 | 49.12 |
2023-11-07 | Tue | 8.16 | 22.14 | 0.00 | 19.35 | 1.22 | 50.86 |
p <- ggplot2::ggplot(recent_dt, aes(x = dv_hh_start_date,
y = sum_kWh_gas, fill = dv_peakPeriod)) +
geom_col(position = "stack") +
scale_fill_viridis_d(name = "Time of day") +
scale_x_date(date_breaks = "1 day", date_labels = "%a %b %d") +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Date",
y = "Gas kWh")
plotly::ggplotly(p)
Figure 5.27: Half hourly gas import (current year, last 14 days)
Compare with temperatures from https://www.accuweather.com/en/gb/framlingham/ip13-9/weather-forecast/330898
Figure 5.28 shows the mean daily kWh import with a smoothed curve.
To do: mark weekends etc
gasCons_dt[, dv_month := lubridate::month(dv_hh_start_date, label = TRUE)]
gasCons_dt[, dv_year := lubridate::year(dv_hh_start_date)]
gasCons_dt[, dv_yday := lubridate::yday(dv_hh_start_date)]
plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_yday, dv_month, dv_year)]
makeDailyPlotByYear(plotDT, yVar = "sum_kWh")
## Scale for alpha is already present.
## Adding another scale for alpha, which will replace the existing scale.
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.28: Daily gas consumption (current year)
plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_hh_start_date, dv_weekend)]
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = sum_kWh,
colour = dv_weekend)) +
geom_point() +
geom_smooth() +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Weekend") +
labs(x = "Date",
y = "Sum kWh per day")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.29: Daily gas consumption (by weekday, current year)
Repeat for mean
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = mean_kWh,
colour = dv_weekend)) +
geom_point() +
geom_smooth() +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Weekend") +
labs(x = "Date",
y = "Mean kWh per day")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.30 shows the mean daily kWh import with a smoothed curve by period of the day.
To do: mark weekends etc
plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_hh_start_date, dv_peakPeriod, dv_weekend)]
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = mean_kWh,
colour = dv_peakPeriod)) +
geom_line() +
geom_smooth() +
#facet_grid(dv_peakPeriod ~ .) +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Peak period") +
labs(x = "Date",
y = "Mean kWh per period")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.30: Daily gas consumption by peak period (current year)
This is much more simple. We can only apply the BEIS 2021 value as there are no time-varying emissions factors for gas.
rmdParams$BEIS_gas_ci <- 0.20297
As before, for the BEIS method we’ll have to use the 2021 emissions factor as the 2022 value is not yet available.
For 2021 this is: 0.20297 Kg CO2e/kWh
gasCons_dt[, KgCO2_beis := consumption_kWh * rmdParams$BEIS_gas_ci]
t <- gasCons_dt[, .(sumkWh = sum(consumption_kWh),
sumKgCO2_beis = sum(KgCO2_beis))]
makeFlexTable(t, cap = "Emissions estimation using gas kWh to date")
sumkWh | sumKgCO2_beis |
|---|---|
21,106 | 4,284 |
Analyse costs using:
The latter are (currently) the same as the at UK price cap: £0.1031 / kWh & £0.2684 ( see Ofgem)
Yes, I know I can extract our exact tariff from the octopus API…
daily_gas <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh, na.rm = TRUE), # beware missing (N/A) may decrease sum
nObs = .N), keyby = .(dv_hh_start_date)]
# extract from pricesDT
# must be an easier way
daily_gas[, kwh_p := ifelse(dv_hh_start_date < lubridate::as_date("2022-10-04"),
pricesDT[fuel == "gas" & component == "kWh" &
dateEnd %like% "2022-10-03", price], # why does it need to be like??
pricesDT[fuel == "gas" & component == "kWh" &
dateStart %like% "2022-10-04", price])]
daily_gas[, sc_p := ifelse(dv_hh_start_date < lubridate::as_date("2022-10-04"),
pricesDT[fuel == "gas" & component == "sc" &
dateEnd %like% "2022-10-03", price],
pricesDT[fuel == "gas" & component == "sc" &
dateStart %like% "2022-10-04", price])]
daily_gas[, cost := ((sum_kWh * kwh_p) + sc_p)]
daily_gas[, month := lubridate::month(dv_hh_start_date, label = TRUE)]
ggplot2::ggplot(daily_gas, aes(x = dv_hh_start_date, y = cost)) +
geom_point(aes(colour = month)) +
geom_smooth() +
geom_vline(xintercept = lubridate::as_date("2022-10-04")) +
labs(y = "Gas daily cost £",
caption = "Tariff change/price cap/EPG 2022 date shown\nSmoothed within month")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Figure 5.31: Daily gas costs
lastWeek <- max(daily_gas$dv_hh_start_date) - 7
makeFlexTable(daily_gas[dv_hh_start_date > lastWeek, .(dv_hh_start_date,
day = lubridate::wday(dv_hh_start_date, label = TRUE),
sum_kWh, nObs, cost)], digits = 2,
cap = "Recent daily gas cost")
dv_hh_start_date | day | sum_kWh | nObs | cost |
|---|---|---|---|---|
2023-11-01 | Wed | 39.48 | 48 | 4.09 |
2023-11-02 | Thu | 68.02 | 47 | 6.87 |
2023-11-03 | Fri | 55.85 | 46 | 5.68 |
2023-11-04 | Sat | 51.55 | 48 | 5.27 |
2023-11-05 | Sun | 50.76 | 47 | 5.19 |
2023-11-06 | Mon | 49.12 | 48 | 5.03 |
2023-11-07 | Tue | 50.86 | 47 | 5.20 |
daily_gas[, month_floor := lubridate::floor_date(dv_hh_start_date, "months")]
monthly_gas <- daily_gas[, .(sum_kWh = sum(sum_kWh),
cost = sum(cost)),
keyby = .(month_floor)]
ggplot2::ggplot(monthly_gas, aes(x = month_floor, y = cost)) +
geom_col() +
geom_vline(xintercept = lubridate::as_date("2022-10-04")) +
labs(y = "Monthly cost £",
x = "Month",
caption = "Tariff change/price cap/EPG 2022 date shown\nBeware incomplete months")
Figure 5.32: Monthly gas costs
message("Projected annual gas total kWh")
## Projected annual gas total kWh
projAnnual_gas_kWh <- mean(daily_gas$sum_kWh)*365
projAnnual_gas_kWh
## [1] 14617.8
message("########")
## ########
message("# Prices to 10th October 2022")
## # Prices to 10th October 2022
# get_price <- function(dt, fuel, component, dateEnd){
# p <- dt[fuel == fuel &
# component == component &
# dateEnd == dateEnd,
# price]
# return(p)
# }
kWh_p <- pricesDT[fuel == "gas" & component == "kWh" &
dateEnd %like% "2022-10-03", price]
sc_p <- pricesDT[fuel == "gas" & component == "sc" &
dateEnd %like% "2022-10-03", price]
message("Projected annual gas cost @ £ ", kWh_p,
" per kWh & standing charge at £ ", sc_p,
" per day")
## Projected annual gas cost @ £ 0.0719 per kWh & standing charge at £ 0.261 per day
projAannualCost <- (projAnnual_gas_kWh*kWh_p)+(365*sc_p) # standing charge
projAannualCost
## [1] 1146.285
message("Mean projected monthly cost: £")
## Mean projected monthly cost: £
projAannualCost/12
## [1] 95.52373
message("########")
## ########
message("# From 4th Oct 2022 - price cap")
## # From 4th Oct 2022 - price cap
kWh_p <- pricesDT[fuel == "gas" & component == "kWh" &
dateStart %like% "2022-10-04", price]
sc_p <- pricesDT[fuel == "gas" & component == "sc" &
dateStart %like% "2022-10-04", price]
message("Projected annual gas cost @ £ ", kWh_p,
" per kWh & standing charge at £ ", sc_p,
" per day")
## Projected annual gas cost @ £ 0.0971 per kWh & standing charge at £ 0.2616 per day
annualCost_gasCapped <- (projAnnual_gas_kWh*kWh_p)+(365*sc_p) # standing charge
annualCost_gasCapped
## [1] 1514.872
message("Projected mean monthly £ under price cap")
## Projected mean monthly £ under price cap
monthlyCost_capped <- annualCost_gasCapped/12
monthlyCost_capped
## [1] 126.2394
message("That's an increase of ", round(100*((annualCost_gasCapped-projAannualCost)/projAannualCost),2), " % points")
## That's an increase of 32.15 % points
Partly in response to: https://twitter.com/heatpolicyrich/status/1603768675223437319
Table 5.22 shows recent usage and costs (see pricing detail above) while Figure ?? shows total cost over time.
daily_gas[, gas_sum_kWh := sum_kWh]
daily_gas[, gas_cost := cost]
daily_elec[, elec_sum_kWh := sum_kWh]
daily_elec[, elec_cost := cost]
daily_gas[, dow := lubridate::wday(dv_hh_start_date, label = TRUE)]
daily_costs <- daily_gas[, .(dv_hh_start_date, dow, gas_sum_kWh, gas_cost, month)][daily_elec[,.(dv_hh_start_date, elec_sum_kWh, elec_cost)]]
daily_costs[, total_kWh := gas_sum_kWh + elec_sum_kWh]
daily_costs[, total_cost := gas_cost + elec_cost]
t <- tail(daily_costs, 10)
makeFlexTable(t, cap = "Latest daily use & cost", digits = 2)
dv_hh_start_date | dow | gas_sum_kWh | gas_cost | month | elec_sum_kWh | elec_cost | total_kWh | total_cost |
|---|---|---|---|---|---|---|---|---|
2023-10-29 | Sun | 34.72 | 3.63 | Oct | 12.99 | 4.60 | 47.70 | 8.23 |
2023-10-30 | Mon | 49.38 | 5.06 | Oct | 13.09 | 4.64 | 62.47 | 9.69 |
2023-10-31 | Tue | 54.83 | 5.59 | Oct | 8.31 | 3.09 | 63.15 | 8.68 |
2023-11-01 | Wed | 39.48 | 4.09 | Nov | 11.91 | 4.25 | 51.38 | 8.35 |
2023-11-02 | Thu | 68.02 | 6.87 | Nov | 9.35 | 3.43 | 77.38 | 10.29 |
2023-11-03 | Fri | 55.85 | 5.68 | Nov | 9.52 | 3.48 | 65.37 | 9.17 |
2023-11-04 | Sat | 51.55 | 5.27 | Nov | 12.37 | 4.40 | 63.92 | 9.67 |
2023-11-05 | Sun | 50.76 | 5.19 | Nov | 7.91 | 2.96 | 58.66 | 8.15 |
2023-11-06 | Mon | 49.12 | 5.03 | Nov | 10.46 | 3.79 | 59.58 | 8.82 |
2023-11-07 | Tue | 50.86 | 5.20 | Nov | 9.34 | 3.42 | 60.20 | 8.62 |
The following two plots are both stacked so the height of the columns shows daily totals. Figure 5.33 shows total daily energy costs.
plotDT <- melt(daily_costs[, .(dv_hh_start_date, gas_cost, elec_cost)],
id.vars = "dv_hh_start_date")
plotDT[, variable := ifelse(variable == "gas_cost",
"Gas £",
"Electricity £")]
plotDT[, variable := factor(variable, levels=c('Electricity £', 'Gas £'))] # re-order so electricity first
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = value, fill = variable)) +
geom_col(position = "stack") +
geom_vline(aes(xintercept = lubridate::as_date("2022-10-04"), colour = "Gas EPG applied")) +
geom_vline(aes(xintercept = lubridate::as_date("2022-11-21"), colour = "Electricity EPG applied")) +
scale_color_discrete(name = "EPG dates") +
scale_fill_discrete(name = "Energy source") +
labs(x = "Date",
y = "Gas + electricity daily cost £",
caption = "Tariff change/price cap/EPG 2022 date shown")
## Warning: Removed 23 rows containing missing values (`position_stack()`).
Figure 5.33: Total daily electricity & gas costs
We’re definitely using too much gas for hot water… (no heating in summer)
Interesting missing gas data
Figure 5.34 shows total daily kWh.
plotDT <- melt(daily_costs[, .(dv_hh_start_date, gas_sum_kWh, elec_sum_kWh)],
id.vars = "dv_hh_start_date")
plotDT[, variable := ifelse(variable == "gas_sum_kWh",
"Gas kWh",
"Electricity kWh")]
plotDT[, variable := factor(variable, levels=c('Electricity kWh', 'Gas kWh'))] # re-order so electricity first
ggplot2::ggplot(plotDT, aes(x = dv_hh_start_date, y = value, fill = variable)) +
geom_col(position = "stack") +
geom_vline(aes(xintercept = lubridate::as_date("2022-10-04"), colour = "Gas EPG applied")) +
geom_vline(aes(xintercept = lubridate::as_date("2022-11-21"), colour = "Electricity EPG applied")) +
scale_color_discrete(name = "EPG dates") +
scale_fill_discrete(name = "Energy source") +
labs(x = "Date",
y = "Gas + electricity daily kWh",
caption = "")
## Warning: Removed 23 rows containing missing values (`position_stack()`).
Figure 5.34: Total daily electricity & gas kWh
ggplot2::ggplot(plotDT[dv_hh_start_date > Sys.Date() - 14], aes(x = dv_hh_start_date, y = value, colour = variable)) +
geom_line() +
scale_colour_discrete(name = "Energy source") +
labs(x = "Date",
y = "Gas + electricity daily kWh",
caption = "")
Figure 5.35: Total daily electricity & gas kWh (last 14 days)
Monthly estimated costs (averaged over 12 months to guide monthly payment value).
These use the latest Octopus flexible tariff from Oct 2022 which are protected by the government’s Energy Price Guarantee:
t <- pricesDT[dateStart > as.Date("2022-10-01")]
t[, dateStart := as.Date(dateStart)]
makeFlexTable(t, cap = "Price cap prices", digits = 3)
fuel | component | dateStart | dateEnd | price | notes |
|---|---|---|---|---|---|
elec_imp | kWh | 2022-11-22 | 2023-06-30 00:00:00 | 0.324 | updated using statement of 5th May 23 - includes EPG |
elec_imp | sc | 2022-11-22 | 2023-06-30 00:00:00 | 0.400 | updated using statement of 5th May 23 |
gas | kWh | 2022-10-04 | 2023-06-30 00:00:00 | 0.097 | updated using statement of 5th May 23 - includes EPG |
gas | sc | 2022-10-04 | 2023-06-30 00:00:00 | 0.262 | updated using statement of 5th May 23 |
elec_imp | kWh | 2023-07-01 | 0.307 | Notice 20/6/23 | |
elec_imp | sc | 2023-07-01 | 0.420 | Notice 20/6/23 claimed this was current £ but it's not! | |
elec_exp | kWh | 2022-11-22 | 0.114 | statement of 5th May 23 | |
gas | kWh | 2023-07-01 | 0.074 | Notice 20/6/23 | |
gas | sc | 2023-07-01 | 0.275 | Notice 20/6/23 |
According to our EPC:
First we’ll compare with hot water…
message("Assume our August-September gas use is just for hot water (almost certainly always true - no gas hob and no heating on)")
## Assume our August-September gas use is just for hot water (almost certainly always true - no gas hob and no heating on)
hw <- gasCons_dt[dv_hh_start_date >= as.Date("2022-08-01") &
dv_hh_start_date < as.Date("2022-10-01"), .(daily_kWh = sum(consumption_kWh),
nObs = .N),
keyby = .(dv_hh_start_date)]
message("Mean daily gas use in this period = ", round(mean(hw$daily_kWh),2), " kWh")
## Mean daily gas use in this period = 31.92 kWh
annualHW <- 365 * mean(hw$daily_kWh)
message("So estimated annual gas for hot water = ", round(annualHW), " kWh (assuming constant hot water use all year round)")
## So estimated annual gas for hot water = 11652 kWh (assuming constant hot water use all year round)
message("That's ", round(annualHW/2351,2), " times the EPC estimate of 2351 kWh... ")
## That's 4.96 times the EPC estimate of 2351 kWh...
And now heating…
message("Assume the estimate of hot water kWh is true")
## Assume the estimate of hot water kWh is true
message("So estimated annual gas for heating = ", round(projAnnual_gas_kWh - annualHW,2), " kWh")
## So estimated annual gas for heating = 2965.63 kWh
message("That's ", round((projAnnual_gas_kWh - annualHW)/6511,2), " times the EPC estimate of 6511 kWh... ")
## That's 0.46 times the EPC estimate of 6511 kWh...
Now switching to overall costs…
message("# Total projected energy")
## # Total projected energy
message(round(projAnnual_gas_kWh + projAnnual_elec_kWh), " kWh")
## 17960 kWh
message("This is ", round(100*((projAnnual_gas_kWh + projAnnual_elec_kWh)/(65*163))), " % of our EPC 'primary energy'")
## This is 170 % of our EPC 'primary energy'
message("#####")
## #####
message("# Total costs")
## # Total costs
message("Monthly total £:", round((annualCost_gasCapped + annualCost_elecCapped)/12))
## Monthly total £:229
annualTotalcapped <- (annualCost_gasCapped + annualCost_elecCapped)
message("Annual total £:", round(annualTotalcapped,2))
## Annual total £:2742.39
If we then subtract the Energy Bill Support Scheme £400 this gives
# subtract £400 Energy Bill Support Scheme
annualTotalcappedAdjusted <- annualTotalcapped - 400
message("Annual adjusted total £:", round(annualTotalcappedAdjusted,2))
## Annual adjusted total £:2342.39
message("Monthly adjusted total £:", round(annualTotalcappedAdjusted/12,2))
## Monthly adjusted total £:195.2
Note that octopus appear to be ‘paying’ this in instalments of £67 by reducing the direct debit value.
NB - this is a simple average across all days/months and takes no account of usage trends (see above). The octopus estimator is smarter than that :-)
Use skmir::skim()to summarise.
skimr::skim(elecCons_dt)
| Name | elecCons_dt |
| Number of rows | 25000 |
| Number of columns | 19 |
| Key | dv_start |
| _______________________ | |
| Column type frequency: | |
| character | 4 |
| Date | 1 |
| difftime | 1 |
| factor | 4 |
| numeric | 8 |
| POSIXct | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| interval_start | 0 | 1 | 20 | 25 | 0 | 25000 | 0 |
| interval_end | 0 | 1 | 20 | 25 | 0 | 25000 | 0 |
| dv_weekend | 0 | 1 | 6 | 8 | 0 | 3 | 0 |
| ba_wd | 0 | 1 | 7 | 7 | 0 | 2 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_hh_start_date | 0 | 1 | 2022-06-05 | 2023-11-07 | 2023-02-20 | 521 |
Variable type: difftime
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_hh_start_hms | 0 | 1 | 0 secs | 84600 secs | 43200 secs | 48 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| dv_month | 0 | 1 | TRUE | 12 | Jul: 2976, Aug: 2976, Oct: 2976, Sep: 2880 |
| dv_peakPeriod | 0 | 1 | FALSE | 5 | Day: 7294, Ear: 6244, Eve: 4168, Lat: 4168 |
| CI_deciles | 0 | 1 | FALSE | 10 | (16: 2574, [27: 2563, (12: 2531, (20: 2528 |
| ba_wday | 0 | 1 | TRUE | 7 | Mon: 3600, Tue: 3600, Sun: 3592, Wed: 3552 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| CARBON_INTENSITY | 0 | 1 | 165.75 | 64.45 | 27.0 | 112.00 | 169.00 | 219.00 | 323.00 | ▅▇▇▇▂ |
| LOW_CARBON_perc | 0 | 1 | 52.78 | 14.94 | 18.9 | 41.00 | 51.00 | 63.70 | 90.20 | ▂▇▇▅▂ |
| RENEWABLE_perc | 0 | 1 | 32.97 | 15.46 | 2.2 | 20.70 | 31.10 | 44.90 | 74.40 | ▃▇▆▅▁ |
| consumption | 0 | 1 | 0.19 | 0.20 | 0.0 | 0.08 | 0.13 | 0.22 | 1.99 | ▇▁▁▁▁ |
| dv_year | 0 | 1 | 2022.60 | 0.49 | 2022.0 | 2022.00 | 2023.00 | 2023.00 | 2023.00 | ▆▁▁▁▇ |
| consumption_kWh | 0 | 1 | 0.19 | 0.20 | 0.0 | 0.08 | 0.13 | 0.22 | 1.99 | ▇▁▁▁▁ |
| dv_yday | 0 | 1 | 198.13 | 94.46 | 1.0 | 131.00 | 208.00 | 273.00 | 365.00 | ▃▃▇▇▅ |
| KgCO2_ngeso | 0 | 1 | 0.03 | 0.04 | 0.0 | 0.01 | 0.02 | 0.04 | 0.42 | ▇▁▁▁▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_start | 0 | 1 | 2022-06-05 04:00:00 | 2023-11-07 23:30:00 | 2023-02-20 13:45:00 | 25000 |
Use skmir::skim()to summarise.
skimr::skim(gasCons_dt)
| Name | gasCons_dt |
| Number of rows | 25000 |
| Number of columns | 14 |
| Key | NULL |
| _______________________ | |
| Column type frequency: | |
| character | 3 |
| Date | 1 |
| difftime | 1 |
| factor | 2 |
| numeric | 6 |
| POSIXct | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| interval_start | 0 | 1 | 20 | 25 | 0 | 25000 | 0 |
| interval_end | 0 | 1 | 20 | 25 | 0 | 25000 | 0 |
| dv_weekend | 0 | 1 | 6 | 8 | 0 | 3 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_hh_start_date | 0 | 1 | 2022-05-07 | 2023-11-07 | 2023-01-30 | 527 |
Variable type: difftime
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_hh_start_hms | 0 | 1 | 0 secs | 84600 secs | 41400 secs | 48 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| dv_month | 0 | 1 | TRUE | 12 | Jul: 2962, Aug: 2955, Jun: 2860, Oct: 2719 |
| dv_peakPeriod | 0 | 1 | FALSE | 5 | Day: 7344, Ear: 6236, Eve: 4186, Lat: 4087 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| consumption | 0 | 1 | 0.07 | 0.18 | 0 | 0 | 0 | 0.01 | 1.21 | ▇▁▁▁▁ |
| dv_year | 0 | 1 | 2022.56 | 0.50 | 2022 | 2022 | 2023 | 2023.00 | 2023.00 | ▆▁▁▁▇ |
| consumption_m3 | 0 | 1 | 0.07 | 0.18 | 0 | 0 | 0 | 0.01 | 1.21 | ▇▁▁▁▁ |
| consumption_kWh | 0 | 1 | 0.84 | 2.03 | 0 | 0 | 0 | 0.12 | 13.79 | ▇▁▁▁▁ |
| dv_yday | 0 | 1 | 199.16 | 90.30 | 1 | 138 | 204 | 269.00 | 365.00 | ▃▅▇▇▅ |
| KgCO2_beis | 0 | 1 | 0.17 | 0.41 | 0 | 0 | 0 | 0.03 | 2.80 | ▇▁▁▁▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_start | 0 | 1 | 2022-05-07 06:30:00 | 2023-11-07 23:00:00 | 2023-01-30 17:15:00 | 25000 |
Figure 6.1 shows the NG-ESO half-hourly carbon intensity over time for the data period as context.
ggplot2::ggplot(elecCons_dt, aes(x = dv_hh_start_date, y = dv_hh_start_hms, fill = CARBON_INTENSITY)) +
geom_tile() +
scale_fill_continuous(name = "Carbon intensity", low = "green", high = "red") +
labs(x = "Date",
y = "Time of day",
caption = "Source: NG-ESO (https://data.nationalgrideso.com/carbon-intensity1/historic-generation-mix)")
Figure 6.1: Half-hourly carbon intensity over time for the data period