DTU02417A4: State Space Model and Kalman Filter

library(lubridate)
library(FKF)
source("KalmanFilterFunction.R")

frameData <- function(){
    dat <- read.csv("Data/data_building.csv")
    vecTimeStamp <- ymd_hms(dat[, 1], "%Y%m%d %H:%M:%S", tz = "UTC")
    numObs <- length(vecTimeStamp) - 1
    dat <- data.frame("series" = seq(1, numObs), "time" = vecTimeStamp[1:numObs], "tempIn" = dat[,2], 
                      "tempAir" = dat[,3], "heat" = dat[,4], "solar" = dat[,5])
    return(dat)
}
dat <- frameData()

 1 failed to parse.

rm(frameData)
numTrain <- length(dat$series) - 3

Question 4.2

plotTrainTestData <- function(vecX_all, vecY_all, numTrain, vecY_test, strNameData, strPosition = "topleft"){
    numObs <- length(vecY_all)
    plot(vecX_all, vecY_all, type = "l", col = "blue", lwd = 1, lty = 1, cex = 0.8, bty = "n", 
         main = paste("Training and Testing Data of", strNameData), 
         xlab = "Series", ylab = strNameData)
    points(seq(numTrain + 1, numObs), vecY_test, col = "blue", lty = 1, pch = 16, cex = 0.8)
    legend(strPosition, inset = .02, legend = c("Training Data", "Testing Data"), col = c("blue", "blue"), 
           pch = c(NaN, 16), lty = c(1, 1), lwd = c(2, 2))
}

plotTrainTestData(dat$series, dat$tempIn, numTrain, tail(dat$tempIn, n = 3), "Indoor Temperature")

plotTrainTestData(dat$series, dat$tempAir, numTrain, tail(dat$tempAir, n = 3), "Outdoor Temperature")

plotTrainTestData(dat$series, dat$solar, numTrain, tail(dat$solar, n = 3), "Solar Irradiation")

plotTrainTestData(dat$series, dat$heat, numTrain, tail(dat$heat, n = 3), "Heating Power")

Question 4.3

matU <- matrix(NA, nrow = 3, ncol = numTrain)
matU[1,] <- dat$tempAir[1:numTrain]
matU[2,] <- dat$heat[1:numTrain]
matU[3,] <- dat$solar[1:numTrain]
matA <- matrix(c(0.775, 0.1, 0.24, 0.9), nrow = 2)       
matB <- matrix(c(0.005, 0, 0.127, 0, 0.335, 0), nrow = 2)      
matC <- matrix(c(1, 0), nrow = 1)
matSigma1_1 <- matrix(c(0.5, 0, 0, 0.5), nrow = 2)
matSigma2_1 <- matrix(0.5)
modKalman_1 <- fkf(a0 = c(25, 25), P0 = diag(c(10, 10), 2), dt = matB %*% matU, ct = 0, Tt = matA, Zt = matC, 
    HHt = matSigma1_1, GGt = matSigma2_1, yt = matrix(dat$tempIn[1:numTrain], nrow = 1), check.input = TRUE)

plotRecons(100, 537, modKalman_1$att[1,], dat$tempIn, modKalman_1$Ptt)

matU.pred <- t(cbind(tail(dat$tempAir, n = 3), tail(dat$heat, n = 3), tail(dat$solar, n = 3)))

predictPlotKalmanFilter(modKalman_1, matSigma1_1, 10, matA, matB, matU.pred)

Question 4.4

resultOptim <- optim(c(25, 25, 10, 10, 0.5, 0.5, 0.5), logLikKalmanFilter, method = "L-BFGS-B", lower = rep(10^-6, 7), 
      upper = c(40, 40, 40, 40, 5, 5, 5))

modKalmanOptim <- calKalmanFilter(resultOptim$par)
matSigma1_optim <- matrix(c(resultOptim$par[5], 0, 0, resultOptim$par[6]), nrow = 2)
predictPlotKalmanFilter(modKalmanOptim, matSigma1_optim, 10, matA, matB, matU.pred)