Domestic hot water forecasting
Peder Bacher
2021-08-17
DHW forecasts
In this vignette we build forecast models for forecasting of Domestic Hot Water (DHW) in different settings.
We will do it in two settings: first forecast the hourly average DHW flow to a multi-family building in Copenhagen, second, forecast the hourly heat load for DHW in single family buildings from Sønderborg.
DHW flow to a multi-family building
The DHW flow was recorded at the supply to a multi-family building in Copenhagen. The data has been resampled to hourly data. It originates from the project “EnergyLab Nordhavn-New Urban Energy Infrastructures and Smart Components” project grant by the Danish Energy Technology Development and Demonstration Program (No. 64015-0055).
Explore and prepare data
First, load the package and data:
# For storing data and cache files
cachedir <- tempdir()
# Load the package and data
library(onlineforecast)
path <- paste0(cachedir,"/DdhwMultiFamily.RDS")
download.file("https://onlineforecasting.org/examples/data/DdhwMultiFamily.RDS", path)
D <- as.data.list(readRDS(path))
# See the variables
str(D)
## List of 2
## $ t : POSIXct[1:2159], format: "2019-01-01 02:00:00" "2019-01-01 03:00:00" ...
## $ dhw: num [1:2159] 412.29 537.45 23.87 227.33 7.99 ...
## - attr(*, "class")= chr "data.list"Plot of all data:
plot_ts(D)
and a zoom of two weeks:
xlim <- c("2019-01-14","2019-01-21")
plot_ts(D, xlim=xlim)
Now generate the time of to be used for a daily curve:
D$tday <- make_tday(D$t, 0:36)
head(D$tday)
## k0 k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0
## 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0 1
## 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0 1 2
## 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0 1 2 3
## 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0 1 2 3 4
## k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36
## 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 3 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 4 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## 5 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 6 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19Make a boxplot conditional on time of day to look for a general pattern:
boxplot(D$dhw ~ D$tday$k0)
Set the score period to have a burn-in of 14 days (needed for a later model):
D$scoreperiod <- in_range(D$t[1]+14*24*3600, D$t)Set up model with single daily curve
Generate a model object:
# New model
model <- forecastmodel$new()
# Clip the output to be inside this range
model$output <- "dhw"
model$outputrange <- c(0,Inf)
# The forgetting factor
model$add_regprm("rls_prm(lambda=0.9)")
model$add_prmbounds(lambda = c(0.9, 0.999, 0.999999))
model$kseq <- 1:36
model$kseqopt <- c(1,18)Add model inputs, first a daily curve with Furrier series basis functions and then an intercept:
model1 <- model$clone_deep()
model1$add_inputs(mu_tday = "fs(tday/24, nharmonics=8)",
mu = "one()")Find the optimal number of harmonics by model stepping:
prm <- list(mu_tday__nharmonics = c(min=2, init=5, max=12))
Lstep1 <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)
##
## ------------------------------------------------------------------------
## Step 1. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=5)
## mu = one()
##
## Current score: 1638.266
## Improvement
## mu_tday__nharmonics=4 "-2.1%"
## mu_tday__nharmonics=6 " 2.7%"
##
## ------------------------------------------------------------------------
## Step 2. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=6)
## mu = one()
##
## Current score: 1594.19
## Improvement
## mu_tday__nharmonics=5 "-2.8%"
## mu_tday__nharmonics=7 " 1.3%"
##
## ------------------------------------------------------------------------
## Step 3. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=7)
## mu = one()
##
## Current score: 1573.916
## Improvement
## mu_tday__nharmonics=6 "-1.3%"
## mu_tday__nharmonics=8 " 0.7%"
##
## ------------------------------------------------------------------------
## Step 4. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=8)
## mu = one()
##
## Current score: 1562.826
## Improvement
## mu_tday__nharmonics=7 "-0.71%"
## mu_tday__nharmonics=9 " 0.14%"
##
## ------------------------------------------------------------------------
## Step 5. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=9)
## mu = one()
##
## Current score: 1560.709
## Improvement
## mu_tday__nharmonics=8 "-0.14%"
## mu_tday__nharmonics=10 " 0.19%"
##
## ------------------------------------------------------------------------
## Step 6. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=10)
## mu = one()
##
## Current score: 1557.724
## Improvement
## mu_tday__nharmonics=9 "-0.19%"
## mu_tday__nharmonics=11 " 0.13%"
##
## ------------------------------------------------------------------------
## Step 7. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=11)
## mu = one()
##
## Current score: 1555.654
## Improvement
## mu_tday__nharmonics=10 "-0.13%"
## mu_tday__nharmonics=12 " 0.19%"
##
## ------------------------------------------------------------------------
## Step 8. Current model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=12)
## mu = one()
##
## Current score: 1552.684
## Improvement
## mu_tday__nharmonics=11 "-0.19%"
##
## ------------------------------------------------------------------------
##
## Final model:
##
## Output: dhw
## Inputs: mu_tday = fs(tday/24, nharmonics=12)
## mu = one()
##
# Keep it
model1 <- Lstep1$final$model
# The trace to see how the level of improvement
getse(Lstep1,"score")
## $step1
## [1] 1638
##
## $step2
## [1] 1594
##
## $step3
## [1] 1574
##
## $step4
## [1] 1563
##
## $step5
## [1] 1561
##
## $step6
## [1] 1558
##
## $step7
## [1] 1556
##
## $final
## [1] 1553Fit the model to all horizons and plot the forecasts:
fit1 <- rls_fit(model1$prm, model1, D)
## ----------------
## lambda
## 1
D$Yhat1 <- fit1$Yhat
plot_ts(D, c("^dhw","^Yhat1"), kseq = c(1,18))
A zoom of a week:
plot_ts(D, c("^dhw|^Yhat1"), xlim, kseq=model$kseqopt)
We see that the last two days, which is weekend, has quite an error, especially the morning peak is later than predicted. So let’s build a model with different diurnal curves between work and weekend days.
First, generate two indicator signals (0 and 1): one indicating workdays and one indicating weekends:
# wday == 0 is Sunday and 6 is Saturday
workday <- as.numeric(aslt(D$t)$wday %in% 1:5)
D$workday <- lagdf(make_input(workday,0:36), 0:-36)
head(D$workday)
## k0 k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
# One for weekend
D$weekend <- 1 - D$workdayplot_ts(D, c("workday","weekend"))
Investigate the residuals of the first model:
Residuals <- D$scoreperiod * residuals(fit1)
par(mfrow=c(3,1))
boxplot(Residuals$h1 ~ D$tday$k0)
boxplot(Residuals$h1 ~ D$tday$k0, subset=D$workday$k0==1)
boxplot(Residuals$h1 ~ D$tday$k0, subset=D$weekend$k0==1)
We can now use these in the model to switch between two curves:
model2 <- model$clone_deep()
model2$add_inputs(mu_tday_workday = "fs(tday/24, nharmonics=10) %**% workday",
mu_tday_weekend = "fs(tday/24, nharmonics=10) %**% weekend",
mu_workday = "workday",
mu_weekend = "weekend")Now optimize and select the two curve model:
# Now optimize
prm <- list(mu_tday_workday__nharmonics = c(min=4, init=12, max=12),
mu_tday_weekend__nharmonics = c(min=4, init=8, max=12))
Lstep2 <- step_optim(model2, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)
##
## ------------------------------------------------------------------------
## Step 1. Current model:
##
## Output: dhw
## Inputs: mu_tday_workday = fs(tday/24, nharmonics=12) %**% workday
## mu_tday_weekend = fs(tday/24, nharmonics=8) %**% weekend
## mu_workday = workday
## mu_weekend = weekend
##
## Current score: 1114.952
## Improvement
## mu_tday_workday__nharmonics=11 "-0.427%"
## mu_tday_weekend__nharmonics=7 " 0.053%"
## mu_tday_weekend__nharmonics=9 "-0.047%"
##
## ------------------------------------------------------------------------
## Step 2. Current model:
##
## Output: dhw
## Inputs: mu_tday_workday = fs(tday/24, nharmonics=12) %**% workday
## mu_tday_weekend = fs(tday/24, nharmonics=7) %**% weekend
## mu_workday = workday
## mu_weekend = weekend
##
## Current score: 1114.357
## Improvement
## mu_tday_workday__nharmonics=11 "-0.428%"
## mu_tday_weekend__nharmonics=6 " 0.042%"
## mu_tday_weekend__nharmonics=8 "-0.053%"
##
## ------------------------------------------------------------------------
## Step 3. Current model:
##
## Output: dhw
## Inputs: mu_tday_workday = fs(tday/24, nharmonics=12) %**% workday
## mu_tday_weekend = fs(tday/24, nharmonics=6) %**% weekend
## mu_workday = workday
## mu_weekend = weekend
##
## Current score: 1113.89
## Improvement
## mu_tday_workday__nharmonics=11 "-0.429%"
## mu_tday_weekend__nharmonics=5 " 0.143%"
## mu_tday_weekend__nharmonics=7 "-0.042%"
##
## ------------------------------------------------------------------------
## Step 4. Current model:
##
## Output: dhw
## Inputs: mu_tday_workday = fs(tday/24, nharmonics=12) %**% workday
## mu_tday_weekend = fs(tday/24, nharmonics=5) %**% weekend
## mu_workday = workday
## mu_weekend = weekend
##
## Current score: 1112.294
## Improvement
## mu_tday_workday__nharmonics=11 "-0.43%"
## mu_tday_weekend__nharmonics=4 " 0.30%"
## mu_tday_weekend__nharmonics=6 "-0.14%"
##
## ------------------------------------------------------------------------
## Step 5. Current model:
##
## Output: dhw
## Inputs: mu_tday_workday = fs(tday/24, nharmonics=12) %**% workday
## mu_tday_weekend = fs(tday/24, nharmonics=4) %**% weekend
## mu_workday = workday
## mu_weekend = weekend
##
## Current score: 1108.929
## Improvement
## mu_tday_workday__nharmonics=11 "-0.43%"
## mu_tday_weekend__nharmonics=5 "-0.30%"
##
## ------------------------------------------------------------------------
##
## Final model:
##
## Output: dhw
## Inputs: mu_tday_workday = fs(tday/24, nharmonics=12) %**% workday
## mu_tday_weekend = fs(tday/24, nharmonics=4) %**% weekend
## mu_workday = workday
## mu_weekend = weekend
##
model2 <- Lstep2$final$model
# Fit to get the forecasts
fit2 <- rls_fit(model2$prm, model2, D)
## ----------------
## lambda
## 0.998
D$Yhat2 <- fit2$YhatPlot the forecasts over a week:
plot_ts(D, c("dhw|Yhat2"), xlim, kseq=model$kseqopt)
Over a month:
plot_ts(D, c("dhw","Yhat2"), xlim=c("2019-02-01","2019-03-01"), kseq=model$kseqopt)
Evaluation
Calculate the score as a function of the horizon and compare the two models:
RMSE <- score(residuals(D[c("Yhat1","Yhat2")], D$dhw), D$scoreperiod)
plot(0, type="n", xlim=range(model$kseq), ylim=range(RMSE), xlab="Horizon k", ylab="RMSE (kW)")
for(i in 1:ncol(RMSE)){
points(model$kseq, RMSE[ ,i], type="b", col=i)
}
legend("right", c("Single daily","Workdays and weekend"), lty=1, col=1:2)
Finally, have a look at the distribution of the residuals conditional on time of day:
Residuals <- D$scoreperiod * residuals(fit2)
par(mfrow=c(3,1))
boxplot(Residuals$h1 ~ D$tday$k0)
boxplot(Residuals$h1 ~ D$tday$k0, subset=D$workday$k0==1)
boxplot(Residuals$h1 ~ D$tday$k0, subset=D$weekend$k0==1)
As a last thing, have a look at how the coefficients change over time in the two curve model. It can be seen how that they are constant in the periods, where their regression inputs are zero, i.e. in the lower weekend coefficients are seen to be constant in the workday periods:
nms <- names(fit2$Lfitval$k1)
D$coefwork <- fit2$Lfitval$k1[grep("workday",nms)]
D$coefweekend <- fit2$Lfitval$k1[grep("weekend",nms)]
plot_ts(D, c("coefwork","coefweekend"), xlim=c("2019-02-01","2019-03-01"), legendcex=0.7, legendspace=15)
Two separate models
Actually, in the case where all variables depends on different periods, it can be computationally lighter to make a model for each period and then combine the forecasts afterwards.
First, a model for the workdays, we simply set the model output data to NA in the period we don’t want to include:
model1$output <- "y"
D$y <- D$dhw
D$y[D$workday$k0 != 1] <- NA
prm <- list(mu_tday__nharmonics = c(min=2, init=5, max=12))
m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model
##
## ------------------------------------------------------------------------
## Step 1. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=5)
## mu = one()
##
## Current score: 1290.359
## Improvement
## mu_tday__nharmonics=4 "-6.8%"
## mu_tday__nharmonics=6 " 7.9%"
##
## ------------------------------------------------------------------------
## Step 2. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=6)
## mu = one()
##
## Current score: 1188.843
## Improvement
## mu_tday__nharmonics=5 "-8.5%"
## mu_tday__nharmonics=7 " 4.6%"
##
## ------------------------------------------------------------------------
## Step 3. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=7)
## mu = one()
##
## Current score: 1133.791
## Improvement
## mu_tday__nharmonics=6 "-4.9%"
## mu_tday__nharmonics=8 " 3.5%"
##
## ------------------------------------------------------------------------
## Step 4. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=8)
## mu = one()
##
## Current score: 1093.894
## Improvement
## mu_tday__nharmonics=7 "-3.6%"
## mu_tday__nharmonics=9 " 1.7%"
##
## ------------------------------------------------------------------------
## Step 5. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=9)
## mu = one()
##
## Current score: 1075.586
## Improvement
## mu_tday__nharmonics=8 "-1.7%"
## mu_tday__nharmonics=10 " 1.5%"
##
## ------------------------------------------------------------------------
## Step 6. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=10)
## mu = one()
##
## Current score: 1059.244
## Improvement
## mu_tday__nharmonics=9 "-1.54%"
## mu_tday__nharmonics=11 " 0.97%"
##
## ------------------------------------------------------------------------
## Step 7. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=11)
## mu = one()
##
## Current score: 1048.99
## Improvement
## mu_tday__nharmonics=10 "-0.98%"
## mu_tday__nharmonics=12 " 0.69%"
##
## ------------------------------------------------------------------------
## Step 8. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=12)
## mu = one()
##
## Current score: 1041.767
## Improvement
## mu_tday__nharmonics=11 "-0.69%"
##
## ------------------------------------------------------------------------
##
## Final model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=12)
## mu = one()
##
fitworkday <- rls_fit(m$prm, m, D)
## ----------------
## lambda
## 0.998Similarly, for the weekends:
D$y <- D$dhw
D$y[D$weekend$k0 != 1] <- NA
m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model
##
## ------------------------------------------------------------------------
## Step 1. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=5)
## mu = one()
##
## Current score: 1265.234
## Improvement
## mu_tday__nharmonics=4 " 0.76%"
## mu_tday__nharmonics=6 "-0.38%"
##
## ------------------------------------------------------------------------
## Step 2. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=4)
## mu = one()
##
## Current score: 1255.642
## Improvement
## mu_tday__nharmonics=3 " 0.72%"
## mu_tday__nharmonics=5 "-0.76%"
##
## ------------------------------------------------------------------------
## Step 3. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=3)
## mu = one()
##
## Current score: 1246.613
## Improvement
## mu_tday__nharmonics=2 "-2.46%"
## mu_tday__nharmonics=4 "-0.72%"
##
## ------------------------------------------------------------------------
##
## Final model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=3)
## mu = one()
##
fit <- rls_fit(m$prm, m, D)
## ----------------
## lambda
## 0.993And combine them and compare them with the two curve model:
# Combine them
i <- D$workday[ ,pst("k",fit$model$kseq)] == 1
i[is.na(i)] <- FALSE
fit$Yhat[i] <- fitworkday$Yhat[i]
# Plot them
D$Yhat <- fit$Yhat
plot_ts(D, c("dhw","^Yhat$","Yhat2"), xlim=c("2019-02-01","2019-03-01"), kseq=model$kseqopt)
So, they are very similar and the score is also very close:
RMSE <- score(residuals(D[c("Yhat","Yhat2")], D$dhw), D$scoreperiod)
plot(0, type="n", xlim=range(model$kseq), ylim=range(RMSE), xlab="Horizon k", ylab="RMSE (kW)")
for(i in 1:ncol(RMSE)){
points(model$kseq, RMSE[ ,i], type="b", col=i)
}
legend("right", c("Two models","One model"), lty=1, col=1:2)
DHW heat load forecasts
Intro
In this section an example is given of using the onlineforecasting package for forecasting domestic hot water heat load.
Data
Data for the forecasting is taken from a data set collected in Sønderborg, Denmark, see Bacher et al. (2013). It comprises heat load measurements for eleven houses use of power for heating domestic hot water. The houses had a tank with hot water connected to district heating. All times are in UTC and the time stamp for average values are set to the end of the time interval.
Load the package:
path <- paste0(cachedir,"/Ddhw11Houses.RDS")
download.file("https://onlineforecasting.org/examples/data/Ddhw11Houses.RDS", path)
D <- as.data.list(readRDS(path))Take a shorter period and plot each individual:
D <- subset(D, c("2010-02-01","2010-07-01"))
plot_ts(D, c(pst("Pdhw",1:11,"$"),"Pdhw$"))
A shorter period:
plot_ts(D, c(pst("Pdhw",1:11,"$"),"Pdhw$"), xlim=c("2010-04-01","2010-04-14"))D$tday <- make_tday(D$t, 0:36)
head(D$tday)
## k0 k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21
## 1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
## 2 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## 3 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0
## 4 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0 1
## 5 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0 1 2
## 6 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 0 1 2 3
## k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36
## 1 23 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 2 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 4 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 5 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## 6 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
D$scoreperiod <- in_range(D$t[1]+14*24*3600, D$t)
# Take the this one, which is made as the mean of all
output <- "Pdhw"# New model
model <- forecastmodel$new()
model$output <- output
# Clip the output to be inside this range
model$outputrange <- c(0,Inf)
# The forgetting factor
model$add_regprm("rls_prm(lambda=0.9)")
model$add_prmbounds(lambda = c(0.9, 0.999, 0.999999))
# The horizons
model$kseq <- 1:36
model$kseqopt <- c(1,18)
# Add inputs
model1 <- model$clone_deep()
model1$add_inputs(mu_tday = "fs(tday/24, nharmonics=8)",
mu = "one()")
# Step to optimize
prm <- list(mu_tday__nharmonics = c(min=2, init=8, max=12))
model1 <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model
##
## ------------------------------------------------------------------------
## Step 1. Current model:
##
## Output: Pdhw
## Inputs: mu_tday = fs(tday/24, nharmonics=8)
## mu = one()
##
## Current score: 4.36956
## Improvement
## mu_tday__nharmonics=7 "-0.016%"
## mu_tday__nharmonics=9 " 0.385%"
##
## ------------------------------------------------------------------------
## Step 2. Current model:
##
## Output: Pdhw
## Inputs: mu_tday = fs(tday/24, nharmonics=9)
## mu = one()
##
## Current score: 4.352746
## Improvement
## mu_tday__nharmonics=8 "-0.386%"
## mu_tday__nharmonics=10 "-0.031%"
##
## ------------------------------------------------------------------------
##
## Final model:
##
## Output: Pdhw
## Inputs: mu_tday = fs(tday/24, nharmonics=9)
## mu = one()
##
#
fit1 <- rls_fit(model1$prm, model1, D)
## ----------------
## lambda
## 0.995
D$Yhat1 <- fit1$Yhat
plot_ts(D, c(pst("^",output),"^Yhat1"), kseq=model$kseqopt)
#
plot_ts(D, c(pst("^",output),"^Yhat1"), xlim=c("2010-04-01","2010-04-14"), kseq=model$kseqopt)
# wday == 0 is Sunday and 6 is Saturday
workday <- as.numeric(aslt(D$t)$wday %in% 1:5)
D$workday <- lagdf(make_input(workday,0:36), 0:-36)
head(D$workday)
## k0 k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
# One for weekend
D$weekend <- 1 - D$workday
## -----------------------------------------------------------------------------
Residuals <- D$scoreperiod * residuals(fit1)
par(mfrow=c(3,1))
boxplot(Residuals$h1 ~ D$tday$k0)
boxplot(Residuals$h1 ~ D$tday$k0, subset=D$workday$k0==1)
boxplot(Residuals$h1 ~ D$tday$k0, subset=D$weekend$k0==1)
D$y <- D[[output]]
model1$output <- "y"
D$y[D$workday$k0 != 1] <- NA
m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model
##
## ------------------------------------------------------------------------
## Step 1. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=8)
## mu = one()
##
## Current score: 4.426977
## Improvement
## mu_tday__nharmonics=7 "-0.051%"
## mu_tday__nharmonics=9 " 0.510%"
##
## ------------------------------------------------------------------------
## Step 2. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=9)
## mu = one()
##
## Current score: 4.404391
## Improvement
## mu_tday__nharmonics=8 "-0.513%"
## mu_tday__nharmonics=10 "-0.085%"
##
## ------------------------------------------------------------------------
##
## Final model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=9)
## mu = one()
##
fitworkday <- rls_fit(m$prm, m, D)
## ----------------
## lambda
## 0.993
# Similarly, for the weekends:
D$y <- D[[output]]
D$y[D$weekend$k0 != 1] <- NA
m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model
##
## ------------------------------------------------------------------------
## Step 1. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=8)
## mu = one()
##
## Current score: 4.163809
## Improvement
## mu_tday__nharmonics=7 " 0.6%"
## mu_tday__nharmonics=9 "-0.3%"
##
## ------------------------------------------------------------------------
## Step 2. Current model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=7)
## mu = one()
##
## Current score: 4.138879
## Improvement
## mu_tday__nharmonics=6 "-0.35%"
## mu_tday__nharmonics=8 "-0.60%"
##
## ------------------------------------------------------------------------
##
## Final model:
##
## Output: y
## Inputs: mu_tday = fs(tday/24, nharmonics=7)
## mu = one()
##
fit <- rls_fit(m$prm, m, D)
## ----------------
## lambda
## 0.988
D$Yhat <- fit$Yhat
# Combine them
i <- D$workday[ ,pst("k",fit$model$kseq)] == 1
i[is.na(i)] <- FALSE
D$Yhat[i] <- fitworkday$Yhat[i]Residuals <- D$scoreperiod * residuals(fit$Yhat, D[[output]])
plot(Residuals$h24)
plot(Residuals$h1)
RMSE <- score(residuals(D[c("Yhat","Yhat1")], D[[output]]), D$scoreperiod)
plot(0, type="n", xlim=range(model$kseq), ylim=range(RMSE), xlab="Horizon k", ylab="RMSE (kW)")
for(i in 1:ncol(RMSE)){
points(model$kseq, RMSE[ ,i], type="b", col=i)
}
legend("right", c("Two models","One model"), lty=1, col=1:2)
Bacher, Peder, Henrik Madsen, Henrik Aalborg Nielsen, and Bengt Perers. 2013. “Short-Term Heat Load Forecasting for Single Family Houses.” Energy and Buildings 65 (0): 101–12. https://doi.org/http://dx.doi.org/10.1016/j.enbuild.2013.04.022.