diff --git a/analysis/GREENGridModel/_loadPowerData.Rmd b/analysis/GREENGridModel/_loadPowerData.Rmd
index 0b0b5d6fe0eff7d485fa378a97d51d5c7447766d..cf0860fa16349d47309840a697bc4beb9dca9c54 100644
--- a/analysis/GREENGridModel/_loadPowerData.Rmd
+++ b/analysis/GREENGridModel/_loadPowerData.Rmd
@@ -60,19 +60,49 @@ setkey(ggHwDT, linkID)
 ```
 
 ```{r hw.summary.overall}
+
+amPeakStart <- hms::as.hms("07:00:00")
+amPeakEnd <- hms::as.hms("09:00:00")
+pmPeakStart <- hms::as.hms("16:00:00")
+pmPeakEnd <- hms::as.hms("20:00:00")
+
+rectAlpha <- 0.1
+vLineAlpha <- 0.4
+vLineCol <- "#0072B2" # http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
+myTextSize <- 3
+
+addPeaks <- function(p){
+  # takes a plot, assumes qHour is hms, adds peak period annotations
+  # assumes you've set yMin & yMax already
+  # breaks with facet_grid, scales = "free" so you have to build seperate plots if you want to do that
+  # there is a complex solution (https://stackoverflow.com/questions/27898651/get-scales-range-of-facets-scales-free) but...
+  p <- p + annotate("rect", xmin = amPeakStart,
+             xmax = amPeakEnd, 
+             ymin = yMin, ymax = yMax, 
+             alpha = rectAlpha, fill = vLineCol)
+  p <- p + annotate("rect", xmin = pmPeakStart,
+             xmax = pmPeakEnd, 
+             ymin = yMin, ymax = yMax, 
+             alpha = rectAlpha, fill = vLineCol)
+  return(p)
+}
+
 plotDT <- ggHhDT[, .(linkID)][ggHwDT][, .(meanW = mean(powerW), # merge on linkID here as we summarise (data.table trick)
                       sdW = sd(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, season)]
 
-ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = season)) +
+p <- ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = season)) +
   geom_step() +
   scale_colour_discrete(name = "Season") +
   #scale_colour_viridis(discrete = TRUE) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption, "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ```{r hw.summary.presenceKids}
@@ -81,16 +111,18 @@ plotDT <- ggHhDT[, .(linkID, presenceKids)][ggHwDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, presenceKids, season)]
 
-ggplot2::ggplot(plotDT[!is.na(presenceKids)], aes(x = qHour, y = meanW/1000, colour = presenceKids)) +
+p <- ggplot2::ggplot(plotDT[!is.na(presenceKids)], aes(x = qHour, y = meanW/1000, colour = presenceKids)) +
   geom_step() +
   scale_colour_discrete(name = "Presence of children") +
   #scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) + 
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption, "\n Peak demand periods shaded")
        )
-
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 # need to count households per category for the actual data plotted
 ggHwIDsDT <- ggHwDT[, .(nObs = .N), keyby = linkID]
 
@@ -109,15 +141,18 @@ plotDT <- ggHhDT[, .(linkID, nBedroomsCat)][ggHwDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, nBedroomsCat, season)]
 
-ggplot2::ggplot(plotDT[!is.na(nBedroomsCat)], aes(x = qHour, y = meanW/1000, colour = nBedroomsCat)) +
+p <- ggplot2::ggplot(plotDT[!is.na(nBedroomsCat)], aes(x = qHour, y = meanW/1000, colour = nBedroomsCat)) +
   geom_step() +
   facet_grid(season ~ .) +
   scale_colour_discrete(name = "N bedrooms") +
   #scale_colour_viridis(discrete = TRUE) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption, "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 
 t <- ggHhDT[ggHwIDsDT][, .('N households' = uniqueN(linkID)),
                        keyby = .('N bedrooms' = nBedroomsCat)]
@@ -133,15 +168,17 @@ plotDT <- ggHhDT[, .(linkID, nPeopleCat)][ggHwDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, nPeopleCat, season)]
 
-ggplot2::ggplot(plotDT[!is.na(nPeopleCat)], aes(x = qHour, y = meanW/1000, colour = nPeopleCat)) +
+p <- ggplot2::ggplot(plotDT[!is.na(nPeopleCat)], aes(x = qHour, y = meanW/1000, colour = nPeopleCat)) +
   geom_step() +
   scale_colour_discrete(name = "N people") +
   #scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
   labs(y = "Mean kW (across households)",
-       x = "Time of day",
-       caption = myCaption
+       x = "Time of day",caption = paste0(myCaption, "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 
 t <- ggHhDT[ggHwIDsDT][, .('N households' = uniqueN(linkID)),
                        keyby = .('N people' = nPeopleCat)]
@@ -157,7 +194,7 @@ plotDT <- ggHwDT[, .(meanW = mean(powerW),
                          nObs = .N), keyby = .(linkID, qHour, season)]
 
 # aggregate plot per hh
-ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
+p <- ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
   geom_step() +
   scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
@@ -165,7 +202,12 @@ ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
        x = "Time of day",
        caption = paste0(myCaption, 
                         "\n Hot Water circuits, n = ", 
-                        uniqueN(ggHwDT$circuit)))
+                        uniqueN(ggHwDT$circuit),
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ```{r hw.season.byLocation}
@@ -178,14 +220,19 @@ plotDT <- ggHhDT[, .(linkID, Location)][ggHwDT][, .(meanW = mean(powerW),
                       nObs = .N), keyby = .(qHour, Location)]
 
 # aggregate plot over all households
-ggplot2::ggplot(plotDT[!is.na(Location)], aes(x = qHour, y = meanW/1000, colour = Location)) +
+p <- ggplot2::ggplot(plotDT[!is.na(Location)], aes(x = qHour, y = meanW/1000, colour = Location)) +
   geom_step() +
   #scale_colour_viridis(discrete = TRUE) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
        caption = paste0(myCaption, 
                         "\n Hot Water circuits: n = ", 
-                        uniqueN(ggHwDT$linkID)))
+                        uniqueN(ggHwDT$linkID),
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ## Heat Pump
@@ -249,14 +296,18 @@ plotDT <- ggHhDT[, .(linkID)][ggHpDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, season)]
 
-ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = season)) +
+p <- ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = season)) +
   geom_step() +
   scale_colour_discrete(name = "Season") +
   #scale_colour_viridis(discrete = TRUE) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption,
+                        "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ```{r hp.summary.charts.presenceKids}
@@ -266,16 +317,21 @@ plotDT <- ggHhDT[, .(linkID, presenceKids)][ggHpDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, presenceKids, season)]
 
-ggplot2::ggplot(plotDT[!is.na(presenceKids)], aes(x = qHour, y = meanW/1000, colour = presenceKids)) +
+p <- ggplot2::ggplot(plotDT[!is.na(presenceKids)], aes(x = qHour, y = meanW/1000, colour = presenceKids)) +
   geom_step() +
   scale_colour_discrete(name = "Presence of children") +
   #scale_colour_viridis(discrete = TRUE) +
-  facet_grid(season ~ .)
+  facet_grid(season ~ .) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
        caption = paste0(myCaption, 
                         "\n Heat pumps, n = ", 
-                        uniqueN(ggHpDT$linkID)))
+                        uniqueN(ggHpDT$linkID),
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
   
 # need to count households per category for the actual data plotted
 ggHpIDsDT <- ggHpDT[, .(nObs = .N), keyby = linkID]
@@ -295,15 +351,19 @@ plotDT <- ggHhDT[, .(linkID, nBedroomsCat)][ggHpDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, nBedroomsCat, season)]
 
-ggplot2::ggplot(plotDT[!is.na(nBedroomsCat)], aes(x = qHour, y = meanW/1000, colour = nBedroomsCat)) +
+p <- ggplot2::ggplot(plotDT[!is.na(nBedroomsCat)], aes(x = qHour, y = meanW/1000, colour = nBedroomsCat)) +
   geom_step() +
   scale_colour_discrete(name = "N bedrooms") +
   #scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption,
+                        "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 
 t <- ggHhDT[ggHpIDsDT][, .('N households' = uniqueN(linkID)),
                        keyby = .('N bedrooms' = nBedroomsCat)]
@@ -319,15 +379,19 @@ plotDT <- ggHhDT[, .(linkID, nPeopleCat)][ggHpDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, nPeopleCat, season)]
 
-ggplot2::ggplot(plotDT[!is.na(nPeopleCat)], aes(x = qHour, y = meanW/1000, colour = nPeopleCat)) +
+p <- ggplot2::ggplot(plotDT[!is.na(nPeopleCat)], aes(x = qHour, y = meanW/1000, colour = nPeopleCat)) +
   geom_step() +
   scale_colour_discrete(name = "N people") +
   #scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption,
+                        "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 
 t <- ggHhDT[ggHpIDsDT][, .('N households' = uniqueN(linkID)),
                        keyby = .('N people' = nPeopleCat)]
@@ -343,7 +407,7 @@ plotDT <- ggHpDT[, .(meanW = mean(powerW),
                          nObs = .N), keyby = .(linkID, qHour, season)]
 
 # aggregate plot per hh
-ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
+p <- ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
   geom_step() +
   scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
@@ -351,7 +415,12 @@ ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
        x = "Time of day",
        caption = paste0(myCaption, 
                         "\n Heat Pump circuits, n = ", 
-                        uniqueN(ggHpDT$circuit)))
+                        uniqueN(ggHpDT$circuit),
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ```{r hp.season.byLocation}
@@ -364,14 +433,19 @@ plotDT <- ggHhDT[, .(linkID, Location)][ggHpDT][, .(meanW = mean(powerW),
                       nObs = .N), keyby = .(qHour, Location)]
 
 # aggregate plot over all households
-ggplot2::ggplot(plotDT[!is.na(Location)], aes(x = qHour, y = meanW/1000, colour = Location)) +
+p <- ggplot2::ggplot(plotDT[!is.na(Location)], aes(x = qHour, y = meanW/1000, colour = Location)) +
   geom_step() +
   #scale_colour_viridis(discrete = TRUE) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
        caption = paste0(myCaption, 
                         "\n Heat Pump circuits: n = ", 
-                        uniqueN(ggHpDT$linkID)))
+                        uniqueN(ggHpDT$linkID),
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ## Lighting
@@ -439,13 +513,18 @@ plotDT <- ggHhDT[, .(linkID)][ggLightingDT][, .(meanW = mean(powerW),
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, season)]
 
-ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = season)) +
+p <- ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = season)) +
   geom_step() +
   scale_colour_discrete(name = "Season") +
   #scale_colour_viridis(discrete = TRUE) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption)
+       caption = paste0(myCaption,
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ```{r l.summary.charts.presenceKids}
@@ -455,14 +534,19 @@ plotDT <- ggHhDT[, .(linkID, presenceKids)][ggLightingDT][, .(meanW = mean(power
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, presenceKids, season)]
 
-ggplot2::ggplot(plotDT[!is.na(presenceKids)], aes(x = qHour, y = meanW/1000, colour = presenceKids)) +
+p <- ggplot2::ggplot(plotDT[!is.na(presenceKids)], aes(x = qHour, y = meanW/1000, colour = presenceKids)) +
   geom_step() +
   scale_colour_discrete(name = "Presence of children") +
   #scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption)
+       caption = paste0(myCaption,
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
   
 # need to count households per category for the actual data plotted
 ggLightingIDsDT <- ggLightingDT[, .(nObs = .N), keyby = linkID]
@@ -482,15 +566,19 @@ plotDT <- ggHhDT[, .(linkID, nBedroomsCat)][ggLightingDT][, .(meanW = mean(power
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, nBedroomsCat, season)]
 
-ggplot2::ggplot(plotDT[!is.na(nBedroomsCat)], aes(x = qHour, y = meanW/1000, colour = nBedroomsCat)) +
+p <- ggplot2::ggplot(plotDT[!is.na(nBedroomsCat)], aes(x = qHour, y = meanW/1000, colour = nBedroomsCat)) +
   geom_step() +
   scale_colour_discrete(name = "N bedrooms") +
   #scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption,
+                        "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 
 t <- ggHhDT[ggLightingIDsDT][, .('N households' = uniqueN(linkID)),
                        keyby = .('N bedrooms' = nBedroomsCat)]
@@ -506,15 +594,19 @@ plotDT <- ggHhDT[, .(linkID, nPeopleCat)][ggLightingDT][, .(meanW = mean(powerW)
                       sumW = sum(powerW),
                       nObs = .N), keyby = .(qHour, nPeopleCat, season)]
 
-ggplot2::ggplot(plotDT[!is.na(nPeopleCat)], aes(x = qHour, y = meanW/1000, colour = nPeopleCat)) +
+p <- ggplot2::ggplot(plotDT[!is.na(nPeopleCat)], aes(x = qHour, y = meanW/1000, colour = nPeopleCat)) +
   geom_step() +
   scale_colour_discrete(name = "N people") +
   #scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
-       caption = myCaption
+       caption = paste0(myCaption,
+                        "\n Peak demand periods shaded")
        )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 
 t <- ggHhDT[ggLightingIDsDT][, .('N households' = uniqueN(linkID)),
                        keyby = .('N people' = nPeopleCat)]
@@ -531,7 +623,7 @@ plotDT <- ggLightingDT[, .(meanW = mean(powerW),
                          nObs = .N), keyby = .(linkID, qHour, season)]
 
 # aggregate plot per hh
-ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
+p <- ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
   geom_step() +
   scale_colour_viridis(discrete = TRUE) +
   facet_grid(season ~ .) +
@@ -539,7 +631,12 @@ ggplot2::ggplot(plotDT, aes(x = qHour, y = meanW/1000, colour = linkID)) +
        x = "Time of day",
        caption = paste0(myCaption, 
                         "\n Lighting circuits, n = ", 
-                        uniqueN(ggLightingDT$circuit)))
+                        uniqueN(ggLightingDT$circuit),
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ```{r l.season.byLocation}
@@ -552,14 +649,19 @@ plotDT <- ggHhDT[, .(linkID, Location)][ggLightingDT][, .(meanW = mean(powerW),
                       nObs = .N), keyby = .(qHour, Location)]
 
 # aggregate plot over all households
-ggplot2::ggplot(plotDT[!is.na(Location)], aes(x = qHour, y = meanW/1000, colour = Location)) +
+p <- ggplot2::ggplot(plotDT[!is.na(Location)], aes(x = qHour, y = meanW/1000, colour = Location)) +
   geom_step() +
   #scale_colour_viridis(discrete = TRUE) +
   labs(y = "Mean kW (across households)",
        x = "Time of day",
        caption = paste0(myCaption, 
                         "\n Lighting circuits: n = ", 
-                        uniqueN(ggLightingDT$linkID)))
+                        uniqueN(ggLightingDT$linkID),
+                        "\n Peak demand periods shaded")
+       )
+yMin <- min(plotDT$meankW)
+yMax <- max(plotDT$meankW)
+addPeaks(p)
 ```
 
 ```{r runToHere}