library(reticulate) library(shiny) library(dplyr) library(ggplot2) library(shinycssloaders) library(TSA) theme_set(theme_bw()) use_virtualenv("../venv/") p <- import("pandas") aggdf <- p$read_pickle("../data/9-clusters.agg.pkl") aggdf$cluster <- factor(aggdf$cluster) sns <- import("seaborn") models <- readRDS('../models/1kmods.rds') clusters = names(models) cbp <- as.character(p$Series(sns$color_palette("colorblind", as.integer(length(clusters)))$as_hex())) names(cbp) <- clusters minhwm <- readRDS("../models/weatherMinharmonic.rds") maxhwm <- readRDS("../models/weatherMaxharmonic.rds") monthchoices <- 1:12 names(monthchoices) <- month.name weekdaychoices = 1:7 names(weekdaychoices) = c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday") ui <- navbarPage("Counties Power demand modelling", collapsible = TRUE, selected = "Prediction", tabPanel("Data and model", plotOutput("dataPlot") %>% withSpinner(type = 5), fluidRow( column(4, selectInput("selDR", "Cluster", choices = clusters, selected = clusters[1]) ), column(4, dateRangeInput("dpDR", "Date range", start = min(as.Date(aggdf$read_time)), end = max(as.Date(aggdf$read_time)), min = min(as.Date(aggdf$read_time)), max = max(as.Date(aggdf$read_time)), format = "dd/mm/yyyy") ), column(4, p(paste("Point represents mean value of cluster, coloured line represents model prediction of", "mean value of cluster.")) ) )), tabPanel("Prediction", plotOutput("predPlot"), fluidRow( column(6, sliderInput("temprange", "Temperature range", min = 0, max = 30, value = c(10, 18), width = "100%", post = " °C", step = 0.5, dragRange = TRUE)), column(2, selectInput("predmon", "Month (approximate)", monthchoices, selected = 1, width = "100%") ), column(2, selectInput("predday", "Day of week", weekdaychoices, selected = 2, width = "100%") ), column(2, numericInput("prl", "Reference line (kwh)", value = 2, width = "100%") ) ), h5("Number of ICPs of each cluster serviced by node:"), fluidRow( column(1, numericInput("c1v", "1:", value = 1, min = 0, step = 1, width = "100%"), offset = 0), column(1, numericInput("c2v", "2:", value = 0, min = 0, step = 1, width = "100%")), column(1, numericInput("c3v", "3:", value = 0, min = 0, step = 1, width = "100%")), column(1, numericInput("c4v", "4:", value = 0, min = 0, step = 1, width = "100%")), column(1, numericInput("c5v", "5:", value = 0, min = 0, step = 1, width = "100%")), column(1, numericInput("c6v", "6:", value = 0, min = 0, step = 1, width = "100%")), column(1, numericInput("c7v", "7:", value = 0, min = 0, step = 1, width = "100%")), column(1, numericInput("c8v", "8:", value = 0, min = 0, step = 1, width = "100%")), column(1, numericInput("c9v", "9:", value = 0, min = 0, step = 1, width = "100%")) ) ) ) server <- function(input, output) { output$dataPlot <- renderPlot({ clus <- input$selDR fildf <- filter(aggdf, cluster %in% clus) %>% mutate(r = resid(models[[clus]]), f = fitted(models[[clus]])) %>% filter(read_time >= as.POSIXct(input$dpDR[1], tz = "UTC"), read_time <= as.POSIXct(input$dpDR[2], tz = "UTC")) yrange = range(c(range(fildf$f), range(fildf$kwh_tot_mean))) dplot <- ggplot(fildf, aes(x = read_time, y = f, color = cluster)) dplot <- dplot + geom_line(size = 2) dplot <- dplot + geom_point(aes(y = kwh_tot_mean), color = "black") dplot <- dplot + coord_cartesian(xlim = as.POSIXct(input$dpDR, tz = "UTC"), ylim = c(yrange[1] - (yrange[2] - yrange[1]) * 0.05, yrange[2] + (yrange[2] - yrange[1]) * 0.05), expand = FALSE) dplot <- dplot + scale_color_manual(values = cbp) dplot <- dplot + theme(legend.position = "none") dplot <- dplot + labs(y = "kwh", x = "Date/Time", caption = sprintf("Cluster: %s", clus)) return(dplot) }) prediction <- function() { numeclus <- c(input$c1v,input$c2v,input$c3v,input$c4v,input$c5v,input$c6v,input$c7v,input$c8v,input$c9v) numeclus[is.na(numeclus)] <- 0 ystart <- floor((as.numeric(input$predmon) - 0.5) / 12 * 365.25 * 48) harm.y <- ts(1:48, frequency = 365.25 * 48, start = c(1, ystart)) %>% harmonic(2) harm.w <- ts(1:48, frequency = 7 * 48, start = c(1, 48 * (-1 + as.numeric(input$predday)))) %>% harmonic(3) harm.d <- ts(1:48, frequency = 48, start = c(1, 1)) %>% harmonic(3) colnames(harm.y) <- sprintf("%s.%s.%s", "year", rep(c("cos", "sin"), each = ncol(harm.y)/2), rep(1:(ncol(harm.y)/2), times = 2)) colnames(harm.w) <- sprintf("%s.%s.%s", "week", rep(c("cos", "sin"), each = ncol(harm.w)/2), rep(1:(ncol(harm.w)/2), times = 2)) colnames(harm.d) <- sprintf("%s.%s.%s", "day", rep(c("cos", "sin"), each = ncol(harm.d)/2), rep(1:(ncol(harm.d)/2), times = 2)) qharm <- harmonic(ts(1:48, frequency = floor(365.25 * 48), start = c(1, ystart)), 2) # harnames <- paste0("yharm", colnames(qharm)) qharm <- as.data.frame(qharm) # names(qharm) <- harnames wmmin <- input$temprange[1] - predict(minhwm, newdata = qharm, type = 'response') wmmax <- input$temprange[2] - predict(maxhwm, newdata = qharm, type = 'response') predinp <- data.frame(resmin = wmmin, resmax = wmmax) predinp <- cbind(predinp, harm.y, harm.w, harm.d) predvec = rep(0, 48) for (c in 1:length(clusters)) { predvec <- predvec + numeclus[c] * predict(models[[c]], newdata = predinp, type = 'response') } return(data.frame(x = 1:48, y = predvec)) } output$predPlot <- renderPlot({ predf <- prediction() predp <- ggplot(predf, aes(x, y)) + geom_line() predp <- predp + scale_y_continuous("Predicted kwh", limits = c(0, NA)) predp <- predp + scale_x_continuous("Time", labels = function(x) { sprintf("%02d:%02d", (x %/% 2) %% 24, x %% 2 * 30) }, breaks = 0:12 * 4) if (!is.na(input$prl)) { predp <- predp + geom_hline(yintercept = input$prl, linetype = "dashed") } predp }) } shinyApp(ui = ui, server = server)