R/utility.R

In msaltieri/termo: Tools for Energy Risk Management & Optimizations

#' Hourly & Quarter-Hourly Calendar Builder
#'
#' This functions builds an hourly or quarter-hourly calendar between two given
#' dates according to the Italian Daylight Saving Time (DST) switching days,
#' i.e. the last Sunday of March (-1 hour) and the last Sunday of October
#' (+1 hour).
#'
#' @param start_date the starting date of the calendar that has to be built.
#' @param end_date the end date of the calendar that has to be built.
#' @param method the method to be used: it could be \code{hourly} (default) or
#' \code{quarter_hourly}.
#' @return \code{calBuilder} returns an object of class \code{data.table}.
#' @author Michele Stefano Altieri, \email{michelestefano.altieri@gmail.com}
#' @keywords calendar dst
#' @examples
#' cal <- calBuilder('2017-01-01', '2017-12-31')
#' cal
#' @import data.table
#' @export
calBuilder <- function(start_date, end_date, method = "hourly") {

    # Questo è per tacitare il check()
    N <- data <- data_ora <- ora <- quarto <- NULL

    # Verifica sugli input
    if (missing(start_date))
        stop("argument `start_date` is missing, with no default")
    if (missing(end_date))
        stop("argument `end_date` is missing, with no default")
    if (!method %in% c("hourly", "quarter_hourly"))
        stop("argument `method` could be choosen between `hourly` and
             `quarter_hourly` only")

    # Converto le date in POSIX
    px_start_date <- as.POSIXct(format(start_date))
    px_end_date <- as.POSIXct(format(as.Date(end_date) + 1)) - 60 * 60
    if (method == "quarter_hourly")
        px_end_date <- px_end_date + 60 * 45

    # Costruisco il calendario orario o quartorario
    cal <- data.table(data_ora = seq(from = px_start_date,
                                     to = px_end_date,
                                     by = ifelse(method == "hourly",
                                                 "1 hour",
                                                 "15 min")))

    # Aggiungo i campi data, ora e quarto (se richiesto)
    sys_tz <- attr(as.POSIXlt(px_start_date), "tzone")[2]
    cal[, data := as.Date(data_ora, tz = sys_tz)]
    cal[, ora := as.integer(1 + hour(data_ora))]
    if (method == "quarter_hourly")
        cal[, quarto := as.integer(1 + 4 * minute(data_ora) / 60)]

    # Correggo il giorno di passaggio da CEST a CET (se presente)
    dstFall <- cal[, .N
                   , by = data][N == 25 * ifelse(method == "hourly", 1, 4),
                                data]
    if (length(dstFall) > 0)
        cal[data %in% dstFall &
                as.POSIXlt(data_ora)$zone == sys_tz,
            ora := ora + 1L]

    # Ritorno il risultato
    return(cal)

}

#' Holidays Builder
#'
#' This function provides a list of every holiday, according with the calendar
#' of the country provided in the input (Italy is the default value but even
#' France or Germany are accepted), returning not-week-end religious and
#' government breaks. Easter and Pentecost days are not returned by the function
#' for they are always on Sunday.
#'
#' @param startDate is the beginning of the calculation period
#' @param endDate is the end of the calculation period
#' @param country is the country for which the holidays are calculated: Italy
#' ('it'), France ('fr') or Germany ('de')
#' @return A vector with the list of the holidays
#' @author Michele Stefano Altieri, \email{michelestefano.altieri@gmail.com}
#' @keywords calendar holidays
#' @examples
#' holidays <- holyBuilder('2016-01-01', '2018-12-31')
#' holidays
#' @importFrom data.table year
#' @export
holyBuilder <- function(startDate, endDate, country = 'it') {

    # Verifica sugli input
    if (missing(startDate))
        stop('argument \'startDate\' is missing, with no default')
    if (missing(endDate))
        stop('argument \'endDate\' is missing, with no default')
    if (!country %in% c('it', 'fr', 'de'))
        stop('input country not recognized')

    # Converto le date nel formato corretto (in caso non lo fossero già)
    startDate <- as.Date(startDate)
    endDate <- as.Date(endDate)

    # Identifico gli anni solari coinvolti
    yearList <- seq(from = year(startDate),
                    to = year(endDate))

    # Definisco la lista dei giorni festivi fissi
    fixedHoly <- switch(country,
                        it = data.frame(mese = c(1, 1, 4, 5, 6, 8, 11, 12, 12, 12),
                                        giorno = c(1, 6, 25, 1, 2, 15, 1, 8, 25, 26)),
                        fr = data.frame(mese = c(1, 5, 5, 7, 8, 11, 11, 12),
                                        giorno = c(1, 1, 8, 14, 15, 1, 11, 25)),
                        de = data.frame(mese = c(1, 1, 5, 8, 10, 11, 12, 12),
                                        giorno = c(1, 6, 1, 15, 3, 1, 25, 26)))

    # Costruisco tutte le festività fisse sugli anni richiesti
    output <- lapply(yearList,
                     function(x) {
                         fixedHolyYear <- data.frame(fixedHoly,
                                                     anno = x)
                         result <- apply(fixedHolyYear, 1, function(y) {
                             return(as.Date(sprintf('%s-%s-%s',
                                                    y['anno'],
                                                    y['mese'],
                                                    y['giorno'])))
                         })
                     })
    output <- as.Date(unlist(output), origin = '1970-01-01')

    # Aggiungo i giorni di Pasquetta
    output <- c(output,as.Date(sapply(yearList,
                                      function(x) {
                                          return(easterBuilder(x) + 1)
                                      }),
                               origin = '1970-01-01'))

    # Se necessario, aggiungo anche il Venerdì Santo
    if (country == 'de')
        output <- c(output,as.Date(sapply(yearList,
                                          function(x) {
                                              return(easterBuilder(x) - 2)
                                          }),
                                   origin = '1970-01-01'))

    # Se necessario, aggiungo anche il giorno dell'Ascensione
    if (country %in% c('fr', 'de'))
        output <- c(output,as.Date(sapply(yearList,
                                          function(x) {
                                              return(easterBuilder(x) + 39)
                                          }),
                                   origin = '1970-01-01'))

    # Se necessario, aggiungo anche la Pentecoste
    if (country %in% c('fr', 'de'))
        output <- c(output,as.Date(sapply(yearList,
                                          function(x) {
                                              return(easterBuilder(x) + 7 * 7 + 1)
                                          }),
                                   origin = '1970-01-01'))

    # Se necessario, aggiungo anche il Corpus Domini
    if (country == 'de')
        output <- c(output,as.Date(sapply(yearList,
                                          function(x) {
                                              return(easterBuilder(x) + 60)
                                          }),
                                   origin = '1970-01-01'))

    # Ordino e filtro il risultato
    output <- sort(output[which(output >= startDate & output <= endDate)])

    # Ritorno il risultato
    return(output)

}

#' Italian Summer Close Period Builder
#'
#' This function provides a list of days that fall in the two central weeks of
#' August. In this period in Italy a large amount of offices, factories and
#' other facilities close, so that both the electric national consumption and
#' the spot prices on the electricity stock market usually suffer a decline.
#'
#' @param startDate is the beginning of the calculation period
#' @param endDate is the end of the calculation period
#' @return A vector with the list of the days in the summer close period
#' @author Michele Stefano Altieri, \email{michelestefano.altieri@gmail.com}
#' @keywords calendar holidays
#' @examples
#' summer <- summerCloseBuilder('2016-01-01', '2018-12-31')
#' summer
#' @importFrom data.table wday year
#' @export
summerCloseBuilder <- function(startDate, endDate) {

    # Verifica sugli input
    if (missing(startDate))
        stop('argument \'startDate\' is missing, with no default')
    if (missing(endDate))
        stop('argument \'endDate\' is missing, with no default')

    # Converto le date nel formato corretto (in caso non lo fossero già)
    startDate <- as.Date(startDate)
    endDate <- as.Date(endDate)

    # Identifico gli anni solari coinvolti
    yearList <- seq(from = year(startDate),
                    to = year(endDate))

    # Costruisco le date nei periodi di chiusura estiva
    output <- lapply(yearList,
                     function(x) {
                         halfAug <- as.Date(ISOdate(x, 8, 15))
                         begin <- halfAug - wday(halfAug - 1) + 1 - 7
                         end <- halfAug - wday(halfAug - 1) + 1 + 6
                         return(seq(begin, end, by = '1 day'))
                     }
    )
    output <- as.Date(unlist(output), origin = '1970-01-01')

    # Ordino e filtro il risultato
    output <- sort(output[which(output >= startDate & output <= endDate)])

    # Ritorno il risultato
    return(output)

}

#' Easter Day Calculator
#'
#' This function computes the Easter Day according with the Gass algorithm.
#' The computation works for years between 1900 and 2099 only!
#'
#' @param year is the year for which the calculation is provided
#' @return the Easter Day
#' @author Michele Stefano Altieri, \email{michelestefano.altieri@gmail.com}
#' @keywords calendar holidays easter
#' @export
easterBuilder <- function(year) {

    a <- year %% 19
    b <- year %% 4
    c <- year %% 7

    d <- (19*a + 24) %% 30
    e <- (2*b + 4*c + 6*d + 5) %% 7

    if ((d + e) < 10) {
        mm <- 3
        dd <- d + e + 22
    } else {
        mm <- 4
        dd <- d + e - 9
    }

    if (mm == 4 & dd == 26) dd <- 19
    if (mm == 4 & dd == 25 & d == 28 & e == 6 & a > 10) dd <- 18

    output <- as.Date(paste(year, mm, dd, sep="-"), "%Y-%m-%d")
    return(output)

}

#' Calendar Time-of-Use Calculator
#'
#' This functions computes the Time-of-Use (TOU) on the provided calendar
#' returning the peak/off-peak international clustering and the one defined by
#' the Italian Authority.
#'
#' @param calendar is a data.table with data_ora, data and ora columns.
#' @param country is the country for which the holidays have to be calculated:
#' Italy (\code{it}) is the default value, but also France (\code{fr}) or
#' Germany (\code{de}) are accepted.
#' @return The function returns the same calendar \code{data.table} with the TOU
#' columns added.
#' @author Michele Stefano Altieri, \email{michelestefano.altieri@gmail.com}
#' @keywords calendar tou peak off-peak
#' @import data.table
#' @export
timeOfUseBuilder <- function(calendar, country = "it") {

    # Questo è per tacitare il check()
    data <- is_holiday <- it2f <- it3f <- ora <- pop <- week_day <- NULL

    # Verifica sull'input
    if (missing(calendar))
        stop("argument `calendar` is missing, with no default")

    # Copio l'input in una nuova struttura
    output <- copy(calendar)

    # Inserisco il giorno della settimana
    output[, week_day := wday(data - 1)]

    # Identifico i giorni di festività
    holidays <- output[, holyBuilder(min(data), max(data), country)]
    output[, is_holiday := data %in% holidays]

    # Calcolo l'aggregazione peak/off-peak
    output[, pop := "OP"]
    output[ora > 8 & ora <= 20 & week_day <= 5,
           pop := "PK"]

    # Calcolo l'aggregazione sulle fasce orarie AEEGSI (3F)
    output[, it3f := "F3"]
    output[(ora %in% 9:19) &
               (week_day %in% 1:5) &
               (!is_holiday),
           it3f := "F1"]
    output[(((ora %in% c(8, 20:23)) &
                 (week_day %in% 1:5)) |
                ((ora %in% 8:23) &
                     (week_day == 6))) &
               (!is_holiday),
           it3f := "F2"]

    # Calcolo l'aggregazione sulle fasce orarie AEEGSI (2F)
    output[, it2f := ifelse(it3f == "F1", "F1", "F23")]

    # Elimino i campi che non servono
    output[, c("week_day", "is_holiday") := list(NULL)]

    # Ritorno il risultato
    return(output)

}

#' Power of a Matrix
#'
#' This function overrides the power operator '^' in order to provide an easy
#' method to perform the power of a square matrix.
#'
#' @param A is the matrix to elevate
#' @param n is the degree of the power
#' @return the n-th power of the A matrix
#' @author Michele Stefano Altieri, \email{michelestefano.altieri@gmail.com}
#' @keywords matrix operations
#' @examples
#' A <- matrix(data = 1:4,
#'             nrow = 2,
#'             ncol = 2)
#' X <- A %^% 2
#' X
#' @export
"%^%" <- function(A, n = 1) {

    # Verifica sugli input
    if (missing(A))
        stop('argument \'A\' is missing, with no default')
    if (round(n) != n)
        stop('argument \'n\' must be integer')

    # Esecuzione del calcolo
    if (n == 0) {

        # Se l'esponente è nullo, ritorno
        # una matrice identità della stessa
        # dimensione dell'input
        return(diag(nrow(A)))

    } else if (n == 1) {

        # Se l'esponente è uno, ritorno
        # la matrice di partenza
        return(A)

    } else {

        # Altrimenti moltiplico n volte
        # la matrice di partenza con se
        # stessa
        B <- A
        for (i in 2:n)
            A <- A %*% B
        return(A)

    }

}