Introducción

El Mars Express Orbiter (MEX) es un satélite de exploración científica de la Agencia Espacial Europea, lanzado en junio del 2003 y colocado en la órbita de Marte. Recolecta información científica sobre la ionósfera y la composición de la superficie del planeta rojo.

El MEX está provisto de siete sensores, entre estos, una cámara estéreo que ha permitido capturar imágenes del planeta. Mucha de la información recolectada por el MEX ha sido crucial para las misiones posteriores de exploración de Marte (Curiosity, Opportunity, etc.).

Figura 1. Imágen con cotas del valle Ancient Atlantis

La nave utiliza energía eléctrica, proveniente de paneles solares (o baterías, durante eclipses) no solo para alimentar los sensores y cámaras, sino también para mantener la temperatura del satélite dentro de rangos óptimos de funcionamiento. El resto de la energía que no se ocupa se puede derivar para realizar las transmisiones de datos.

El objetivo del proyecto consiste en hacer predicción sobre el consumo eléctrico (double, mA) de los 33 circuitos eléctricos de la nave. Al tratarse de una predicción sobre valores discretos, se utilizará un modelo de regresión.

La forma de evaluar la predicción consiste en enviar un archivo de predicción, con promedios por cada hora, para el cuarto año y los 33 circuitos. La precisión de la medición se calculará mediante el Root Mean Square Error (RMSE) que se calcula de la siguiente forma:

\(\epsilon = \sqrt[]{\frac{1}{NM}\sum{(c_{ij}-r_{ij})^{2}}}\)

Para la visualización y la implementación del algoritmo se utilizó R en RStudio.

Debido al volúmen de datos, se trabajará únicamente con los datos del primer año.

Descripción de los datos

Los datasets se obtuvieron desde la página del concurso: https://kelvins.esa.int/mars-express-power-challenge/data/. Están divididos en 4 años marcianos (1 año marciano = 687 días terrestres). Los tres primeros años (2008.08 - 2010.07, 2010.07 - 2012.05, 2012.05 - 2014.04) consisten en datos de entrenamiento. El cuarto año (2014.04 - 2016.03) tiene valores para pruebas.

Cada dataset está a su vez separado según el tipo (o contexto) de información que contienen. Así, existen cinco tipos de contextos de datos: Solar Aspect Angles (SAAF), Detailed Mission Operations Plan (DMOP), Flight Dynamics Timeline (FTL), Other Events(EVTF), Long Term Data such as Sun-Mars Distance (LTDATA). Existe así también un sexto archivo, que contiene los valores de los consumos eléctricos. El cuarto año no tiene este archivo, el cual se debe predecir.

SAAF: Aspectos solares

Los ángulos entregan información sobre la incidencia de los rayos del sol sobre los paneles solares. Los ángulos se miden con respecto a la línea Sol-MEX.

DMOP: Detalles de la planificación de operación

Los archivos DMOP entregan información sobre los comandos que han sido activados en cada sistema. Cada comando incluye el nombre del sistema y el nombre del comando. Debido a la cantidad de comandos, estos no se explican, sin embargo cada uno tiene distintos efectos sobre la temperatura en los distintos subsistemas, y por esto afectando el sistema termal. En estos comandos se encuentran el encendio y apagado del sistema de telecomunicaciones y de los instrumentos científicos.

FTL: Eventos de la trayectoria de la nave

Los eventos listados en estos datos pueden impactar el comportamiento del satélite ya que pueden afectar los ángulos de incidencia, así como también tener relación con el encendido o apagado de algunos instrumentos.

EVTF: Otros eventos

Este set de datos contiene varios eventos. Contiene información sobre eventos del FTL y la complementa con información sobre otros eventos, como por ejemplo eclipses y sus fases (Penumbra y umbra).

LTDATA: Información de periodos extendidos

Incluye información sobre la distancia entre el Sol y Marte, distancia a la Tierra, duración de los eclipses, etc.

POWER: Consumo energético

En estos sets se incluye información sobre la corriente medida en los 33 circuitos eléctricos de la nave. Cada medición se realiza cada 30 o 60 segundos.

Preprocesamiento y Visualización

Datasets

evtf1 = read.csv("Data\\train_set\\context--2008-08-22_2010-07-10--evtf.csv")
ftl1 = read.csv("Data\\train_set\\context--2008-08-22_2010-07-10--ftl.csv")
ltdata1 = read.csv(file="Data\\train_set\\context--2008-08-22_2010-07-10--ltdata.csv")
saaf1 = read.csv(file="Data\\train_set\\context--2008-08-22_2010-07-10--saaf.csv")
dmop1 = read.csv(file="Data\\train_set\\context--2008-08-22_2010-07-10--dmop.csv")
power1 = read.csv(file="Data\\train_set\\power--2008-08-22_2010-07-10.csv")

A continuación se presentan 5 valores de algunas columnas de cada datasets para familiarizar al lector con los datos:

power1[1:5,1:6]
##          ut_ms  NPWD2372   NPWD2401 NPWD2402  NPWD2451    NPWD2471
## 1 1.219363e+12 0.0018210 0.00147389 0.172173 0.0050155 0.000714286
## 2 1.219363e+12 0.0021852 0.00176867 0.177440 0.0060186 0.000714286
## 3 1.219363e+12 0.0021852 0.00147389 0.172173 0.0050155 0.000571429
## 4 1.219363e+12 0.0018210 0.00147389 0.172173 0.0050155 0.000714286
## 5 1.219363e+12 0.0018210 0.00176867 0.177440 0.0050155 0.000714286
dmop1[1:5,]
##          ut_ms subsystem
## 1 1.219366e+12  AXXX380A
## 2 1.219367e+12  ASEQ4200
## 3 1.219367e+12 ATTTF301E
## 4 1.219369e+12 ATTTF310A
## 5 1.219369e+12  APSF01A2
ltdata1[1:5,1:4]
##          ut_ms sunmars_km earthmars_km sunmarsearthangle_deg
## 1 1.219363e+12  241938908    355756044              19.56508
## 2 1.219450e+12  241800160    356303701              19.39007
## 3 1.219536e+12  241660299    356843151              19.21473
## 4 1.219622e+12  241519334    357374355              19.03905
## 5 1.219709e+12  241377272    357897266              18.86303
ftl1[1:5,1:4]
##         utb_ms       ute_ms  type flagcomms
## 1 1.219363e+12 1.219365e+12 EARTH     FALSE
## 2 1.219365e+12 1.219370e+12 EARTH      TRUE
## 3 1.219370e+12 1.219370e+12 EARTH     FALSE
## 4 1.219370e+12 1.219370e+12  SLEW     FALSE
## 5 1.219370e+12 1.219373e+12 NADIR     FALSE
evtf1[1:5,]
##          ut_ms        description
## 1 1.219370e+12         MRB_AOS_00
## 2 1.219370e+12    1200_KM_DESCEND
## 3 1.219371e+12         MRB_AOS_10
## 4 1.219371e+12     800_KM_DESCEND
## 5 1.219371e+12 MAR_PENUMBRA_START
saaf1[1:5,1:4]
##          ut_ms   sa    sx    sy
## 1 1.219363e+12 0.32 14.55 90.32
## 2 1.219363e+12 0.34 14.56 90.34
## 3 1.219363e+12 0.34 14.56 90.34
## 4 1.219363e+12 0.34 14.56 90.34
## 5 1.219363e+12 0.34 14.56 90.34

Escala temporal

Previo a cualquier operación con los datos, se procede a convertir los valores del tiempo.

Originalmente están en milisegundos:

saaf1[1:5,]
##          ut_ms   sa    sx    sy
## 1 1.219363e+12 0.32 14.55 90.32
## 2 1.219363e+12 0.34 14.56 90.34
## 3 1.219363e+12 0.34 14.56 90.34
## 4 1.219363e+12 0.34 14.56 90.34
## 5 1.219363e+12 0.34 14.56 90.34

Para hacer la conversión (https://en.wikipedia.org/wiki/POSIX):

saaf1$ut_ms = as.POSIXct((((saaf1['ut_ms'])/1000)[,]), origin="1970-01-01")

El resultado, en tiempo:

saaf1[1:5,]
##                 ut_ms   sa    sx    sy
## 1 2008-08-21 20:00:13 0.32 14.55 90.32
## 2 2008-08-21 20:00:35 0.34 14.56 90.34
## 3 2008-08-21 20:01:35 0.34 14.56 90.34
## 4 2008-08-21 20:02:35 0.34 14.56 90.34
## 5 2008-08-21 20:03:35 0.34 14.56 90.34

Lo mismo aplica para los demás sets:

ltdata1$ut_ms = as.POSIXct((((ltdata1['ut_ms'])/1000)[,]), origin="1970-01-01")
dmop1$ut_ms = as.POSIXct((((dmop1['ut_ms'])/1000)[,]), origin="1970-01-01")
ftl1$utb_ms = as.POSIXct((((ftl1['utb_ms'])/1000)[,]), origin="1970-01-01")
ftl1$ute_ms = as.POSIXct((((ftl1['ute_ms'])/1000)[,]), origin="1970-01-01")
evtf1$ut_ms = as.POSIXct((((evtf1['ut_ms'])/1000)[,]), origin="1970-01-01")
power1$ut_ms = as.POSIXct((((power1['ut_ms'])/1000)[,]), origin="1970-01-01")

Para varios datasets, se desea comparar el comportamiento de sus variables con respecto al consumo de potencia, por lo que se introduce el siguiente vector:

globalMS = power1$ut_ms

Visualización de POWER: suma + anual

La suma de los consumos para el primer año está dado por la suma de todos los consumos individuales, de la siguiente manera:

sum <- rowSums(power1[,c(2:34)])
ut_ms<-power1$ut_ms
power1Sum <- data.frame(ut_ms,sum)
ggplot(power1Sum,aes(x=ut_ms, y=sum)) + geom_line()

El caso de los consumos individuales se verá más adelante, ya que a esta escala no se puede distinguir.

Visualización de SAAF: sa + anual

Si bien la incidencia solar se compone de las 3 coordenadas (sx,sy,sz) la gráfica total sirve para realizar una visualización y comparación con el consumo de potencia.

 ggplot(saaf1,aes(x=ut_ms, y=sa)) + geom_line()

A partir del gráfico anterior, se puede notar cierta relación directa entre la incidencia solar y el consumo.

Visualización de EVTFL: umbra y penumbra + anual

La umbra y la penumbra son situaciones de oscuridad (parcial para la penumbra). https://es.wikipedia.org/wiki/Penumbra. En estos casos, la intensidad solar disminuirá, por lo que afectará la predicción.

El dataset contiene mucha información de muchos tipos. Particularmente interesan los eventos de umbra y penumbra. Se pueden filtrar aquellos que contienen “UMBRA”:

evtf1Filter = evtf1[grepl("UMBRA", evtf1$description),]
evtf1Filter[1:5,]
##                 ut_ms        description
## 1 2008-08-21 22:08:22 MAR_PENUMBRA_START
## 2 2008-08-21 22:09:46    MAR_UMBRA_START
## 3 2008-08-21 22:11:15      MAR_UMBRA_END
## 4 2008-08-21 22:12:32   MAR_PENUMBRA_END
## 5 2008-08-22 04:59:26 MAR_PENUMBRA_START

Para poder compararlo con el dataset de consumos (POWER), es necesario hacer match sobre las valores temporales:

#Match con la escala de tiempo global:
evtf1Scale = data.frame("ut_ms" = globalMS)
evtf1Filter$ut_ms = globalMS[sapply(evtf1Filter$ut_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]

Para crear rangos con este tipo de eventos, es necesario separar cada columna, según la ocurrencia:

#Transponer:
evtf1Groups = dcast(evtf1Filter, ut_ms ~ description, length)

Posteriormente, hacer un merge con la escala completa de tiempos, rellenar los NA por 0:

#Completar los valores faltantes:
evtf1Scale=merge(x=evtf1Scale, y=evtf1Groups, by="ut_ms", all.x=TRUE)
evtf1Scale[is.na(evtf1Scale)] = 0

Luego se pueden crear pares de variables donde empieza y termina un evento:

#Pares de Penumbra/Umbra
ePairs = data.frame(rbind(c("DEI_PENUMBRA_START","DEI_PENUMBRA_END"), 
                          c("MAR_PENUMBRA_START", "MAR_PENUMBRA_END"),
                          c("MAR_UMBRA_START", "MAR_UMBRA_END"),
                          c("PHO_PENUMBRA_START", "PHO_PENUMBRA_END"),
                          c("PHO_UMBRA_START", "PHO_UMBRA_END")))

Y con estos, crear rangos:

evtf1Ranges = data.frame("ut_ms" = evtf1Scale$ut_ms)

#Iterar sobre los pares:
for (r in 1:nrow(ePairs)) {
  eRange = evtf1Scale[, c(1,match(as.character(ePairs[r,1]),colnames(evtf1Scale)),match(as.character(ePairs[r,2]),colnames(evtf1Scale)))]
  names(eRange)[3] = "end"
  names(eRange)[2] = "start"
  eRange$duration = with(eRange, ave(start-end, 1, FUN = cumsum))
  eRange$duration[eRange$end == 1] = 1
  names(eRange)[4] = paste(as.character(ePairs[r,1]), as.character(ePairs[r,2]), sep = " - ")
  evtf1Ranges = cbind(evtf1Ranges, eRange[4])
}

#Juntar los pares:
evtf1RangesMelt <- melt(evtf1Ranges, id='ut_ms')
ggplot(evtf1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

La gráfica calza parcialmente en rangos con las variaciones del consumo energético, por lo que será considerado tentativamente como una variable influyente.

Visualización de FTL

En el caso de FTL, existen dos tiempos, un tiempo de inicio y un tiempo de término. Se desea construir un dataframe que tenga una sola columna de tiempo y ventanas, al igual que en el caso de la EVTL.

En primera instancia, se visualizará la ocurrencia de las combinaciones de comandos:

ftl1Sys = ddply(ftl1, .(type, flagcomms), nrow)
ftl1Sys = ftl1Sys[ order(-ftl1Sys[,3], ftl1Sys[,1]), ]
  ftl1Sys[1:10,]
##             type flagcomms   V1
## 1           SLEW     FALSE 7267
## 2          EARTH     FALSE 4461
## 3          EARTH      TRUE 2239
## 4       INERTIAL     FALSE 1810
## 5          NADIR     FALSE  929
## 6    MAINTENANCE     FALSE  641
## 7  RADIO_SCIENCE      TRUE  548
## 8         WARMUP     FALSE  474
## 9   ACROSS_TRACK     FALSE  367
## 10        D4PNPO     FALSE  327

De la misma manera, se buscarán aquellos comandos donde las ventanas entre begin y end sean mayores:

ftl1diff = ftl1
ftl1diff$diff = ftl1diff$ute_ms-ftl1diff$utb_ms
ftl1diff$utb_ms = NULL
ftl1diff$ute_ms = NULL

ftl1diffSum = aggregate(ftl1diff$diff, by=list(type=ftl1diff$type, flagcomms=ftl1diff$flagcomms), FUN=sum)
ftl1diffSum = ftl1diffSum[ order(-ftl1diffSum[,3], ftl1diffSum[,1]), ]
  ftl1diffSum[1:10,]
##             type flagcomms        x
## 1          EARTH      TRUE 20186660
## 2          EARTH     FALSE 12578080
## 3           SLEW     FALSE  8889846
## 4       INERTIAL     FALSE  3765371
## 5  RADIO_SCIENCE      TRUE  3509046
## 6    MAINTENANCE     FALSE  3358336
## 7         WARMUP     FALSE  3236773
## 8          NADIR     FALSE  3053966
## 9   ACROSS_TRACK     FALSE  1237024
## 10        WARMUP      TRUE  1184929

De lo anterior se puede apreciar que los 3 type y flagcomms más influyentes son “EARTH,TRUE”, “EARTH,FALSE”, “SLEW,FALSE”. Se filtrará el set para trabajar únicamente con estos valores.

ftl1Filt = ftl1[with(ftl1, grepl("FALSE", ftl1$flagcomms)&grepl("SLEW", ftl1$type)),]
ftl1Filt[1:5,]
##                utb_ms              ute_ms  type flagcomms
## 1 2008-08-21 20:00:13 2008-08-21 20:38:14 EARTH     FALSE
## 2 2008-08-21 21:45:55 2008-08-21 21:46:59 EARTH     FALSE
## 3 2008-08-21 23:09:23 2008-08-21 23:43:53 EARTH     FALSE
## 4 2008-08-22 00:59:04 2008-08-22 02:30:04 EARTH     FALSE
## 5 2008-08-22 02:30:04 2008-08-22 02:31:04 EARTH     FALSE

Para la creación de los rangos:

#Match de los tiempos con el global:
ftl1Filt$utb_ms = globalMS[sapply(ftl1Filt$utb_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]
ftl1Filt$ute_ms = globalMS[sapply(ftl1Filt$ute_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]

#Crear un vector vacío
ftl1Scale = data.frame("ut_ms" = globalMS)
ftl1Scale$type = ""
ftl1Scale$flagcomms = ""
ftl1Scale$start = 0
ftl1Scale$end = 0
#Iterar sobre la colección:
for (i in 1:nrow(ftl1Filt)) {
#for (i in 1:1000) {  
  x = findInterval(ftl1Filt$utb_ms[i], globalMS, all.inside = TRUE, rightmost.closed = TRUE)
  y = findInterval(ftl1Filt$ute_ms[i], globalMS, all.inside = TRUE, rightmost.closed = TRUE)

  if(ftl1Scale[x,]$start == 1 || ftl1Scale[x,]$end == 1 || ftl1Scale[y,]$start == 1 || ftl1Scale[x,]$end == 1){ }
  else{
    ftl1Scale[x,]$type = as.character(ftl1Filt$type[i])
    ftl1Scale[x,]$flagcomms = as.character(ftl1Filt$flagcomms[i])
    ftl1Scale[x,]$start = 1

    ftl1Scale[y,]$type = as.character(ftl1Filt$type[i])
    ftl1Scale[y,]$flagcomms = as.character(ftl1Filt$flagcomms[i])
    ftl1Scale[y,]$end = 1
  }
}

Gráficamente:

ftl1Ranges = data.frame("ut_ms" = ftl1Scale$ut_ms)

ftl1EF = read.csv("Web\\Data\\ftl1Scale-EarthFalse.csv")
fRange = ftl1EF[, c("ut_ms","start","end")]
fRange$duration = with(fRange, ave(start-end, 1, FUN = cumsum))
fRange$duration[fRange$end == 1] = 1
names(fRange)[4] = "EarthFalse"
ftl1Ranges = cbind(ftl1Ranges, fRange[4])

ftl1ET = read.csv("Web\\Data\\ftl1Scale-EarthTrue.csv")
fRange = ftl1ET[, c("ut_ms","start","end")]
fRange$duration = with(fRange, ave(start-end, 1, FUN = cumsum))
fRange$duration[fRange$end == 1] = 1
names(fRange)[4] = "EarthTrue"
ftl1Ranges = cbind(ftl1Ranges, fRange[4])

ftl1RangesMelt <- melt(ftl1Ranges, id='ut_ms')
ggplot(ftl1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

En una escala más cercana:

Visualización de DMOP y LTDATA1: anual

Los comandos del grupo DMOP no se encuentran explicados en la descripción del proyecto. Se puede observar que corresponden a strings.

dmop1[1:5,]
##          ut_ms subsystem
## 1 1.219366e+12  AXXX380A
## 2 1.219367e+12  ASEQ4200
## 3 1.219367e+12 ATTTF301E
## 4 1.219369e+12 ATTTF310A
## 5 1.219369e+12  APSF01A2

La hipótesis es que estos tienen un grupo (group), que corresponde a las 4 primeras letras, un estado (status), que corresponde a los siguientes dos bytes y una cola (end), que corresponde a los dos últimos bytes. Lo que si se indica en el enunciado del problema, es que algunos de estos comandos funcionan como instrucciones ON/OFF a algún circuito de la nave, por lo que su estado podría afectar la predicción.

Luego de examinar los comandos, se ha podido determinar que los comandos importantes son aquellos que empiezan con la letra A, por lo que se procede a eliminar los otros:

#Eliminar los DMOP que no empiezan con "A":
dmop1 = dmop1[grepl("^A", dmop1$subsystem),]

En el siguiente bloque se ordena el código para generar los vectores correspondiendes:

##Generar la columna grupo, según el substring(0,4) del subsystem:
dmop1$group = sapply(dmop1$subsystem, substring, 0, 4)
dmop1$end = sapply(dmop1$subsystem, str_sub, -2,-1)
dmop1$status = sapply(dmop1$subsystem, function(x) paste(sprintf("%04d", strtoi(intToBin(str_sub(x,-4,-4)))),sprintf("%04d", strtoi(intToBin(str_sub(x,-3,-3)))),sep = ""))
dmop1[1:5,]
##                 ut_ms subsystem group end  status
## 1 2008-08-21 20:00:11  AXXX301A  AXXX  1A  110000
## 2 2008-08-21 20:28:29 AAAAF20C1  AAAA  C1  100000
## 3 2008-08-21 20:28:34 AAAAF57A1  AAAA  A1 1010111
## 4 2008-08-21 20:28:39 AAAAF23G1  AAAA  G1  100011
## 5 2008-08-21 20:28:44 AAAAF60A1  AAAA  A1 1100000

Al agrupar los comandos según subsystem, se puede tener una idea de los comandos más influyentes:

##Agregar la cuenta de los elementos:
dmopSys = ddply(dmop1, .(subsystem, group, end, status), nrow)
names(dmopSys)[5] = "Count"
dmopSys = dmopSys[ order(-dmopSys[,5], dmopSys[,1]), ]
dmopSys[1:10,]
##    subsystem group end  status Count
## 1  AAAAF56A1  AAAA  A1 1010110 12054
## 2   APSF28A1  APSF  A1  101000 10582
## 3   AVVV02A0  AVVV  A0      10  8572
## 4   AVVV03B0  AVVV  B0      11  8273
## 5  AAAAF20E1  AAAA  E1  100000  7657
## 6  AAAAF20C1  AAAA  C1  100000  5579
## 7   APSF38A1  APSF  A1  111000  5289
## 8  AAAAF23G1  AAAA  G1  100011  3843
## 9  AAAAF57A1  AAAA  A1 1010111  3279
## 10 AAAAF60A1  AAAA  A1 1100000  3273

Luego de examinar los pares, según la cantidad de apariciones y consideando que el grupo (group) y el final (end) podrían ser candidatos a parejas. De esta forma se crea el siguiente vector con las parejas más influyentes:

#Generar un listado de pares:
dPairs = data.frame(rbind(c("ATTTF310A","ATTTF310B"), 
                          c("AAAAF84A1", "AAAAF85A1"),
                          c("AAAAF57A1", "AAAAF60A1"),
                          c("AAAAF59A1", "AAAAF93A1"),
                          c("ASSSF06A0", "ASSSF01A0"),
                          c("AMMMF19A0", "AMMMF18A0"),
                          c("AMMMF24A0", "AMMMF23A0"),
                          c("AAAAF20D1", "AAAAF18D1"),
                          c("AXXX305B", "AXXX301B"),
                          c("AMMMF20A0", "AMMMF11A0")))
names(dPairs)[1] = "ON"
names(dPairs)[2] = "OFF"

Al filtrar el dataset para aquellos elementos dentro de los pares considerados, se tiene el siguiente frame:

#Filtrar aquellos que estén entre los pares:
dmopFilt = dmop1[which(dmop1$subsystem %in% dPairs$ON | dmop1$subsystem %in%  dPairs$OFF),]
dmopFilt$ut_ms = power1$ut_ms[sapply(dmopFilt$ut_ms, function(x) findInterval(x, power1$ut_ms, all.inside = TRUE, rightmost.closed = TRUE))]
dmopFilt[1:10,]
##                  ut_ms subsystem group end  status
## 1  2008-08-21 20:28:34 AAAAF57A1  AAAA  A1 1010111
## 2  2008-08-21 20:28:44 AAAAF60A1  AAAA  A1 1100000
## 3  2008-08-21 21:37:41 ATTTF310A  ATTT  0A  110001
## 4  2008-08-21 21:53:46 AMMMF11A0  AMMM  A0   10001
## 5  2008-08-21 21:54:01 AMMMF24A0  AMMM  A0  100100
## 6  2008-08-21 21:54:06 AMMMF23A0  AMMM  A0  100011
## 7  2008-08-21 21:54:11 AMMMF20A0  AMMM  A0  100000
## 8  2008-08-21 21:56:36 AMMMF18A0  AMMM  A0   11000
## 9  2008-08-21 21:56:41 AMMMF19A0  AMMM  A0   11001
## 10 2008-08-21 22:41:51 ATTTF310B  ATTT  0B  110001

Al transponer este frame, para tener una relación directa de los tiempos con los comandos, se puede construir el siguiente vector:

#Transponer los filtrados:
dmopGroups = dcast(dmopFilt, ut_ms ~ subsystem, length)
dmop1Scale = data.frame("ut_ms" = globalMS)
dmop1Scale=merge(x=dmop1Scale, y=dmopGroups, by="ut_ms", all.x=TRUE)
dmop1Scale[is.na(dmop1Scale)] = 0
dmopGroups[1:10,c(1,4,6,10,11,12)]
dmopGroups[1:10,]
##                  ut_ms subsystem group end  status
## 1  2008-08-21 20:28:34 AAAAF57A1  AAAA  A1 1010111
## 2  2008-08-21 20:28:44 AAAAF60A1  AAAA  A1 1100000
## 3  2008-08-21 21:37:41 ATTTF310A  ATTT  0A  110001
## 4  2008-08-21 21:53:46 AMMMF11A0  AMMM  A0   10001
## 5  2008-08-21 21:54:01 AMMMF24A0  AMMM  A0  100100
## 6  2008-08-21 21:54:06 AMMMF23A0  AMMM  A0  100011
## 7  2008-08-21 21:54:11 AMMMF20A0  AMMM  A0  100000
## 8  2008-08-21 21:56:36 AMMMF18A0  AMMM  A0   11000
## 9  2008-08-21 21:56:41 AMMMF19A0  AMMM  A0   11001
## 10 2008-08-21 22:41:51 ATTTF310B  ATTT  0B  110001

Al realizar el procedimiento anterior para completar los valores faltantes y graficar, se obtiene el siguiente resultado:

#Crear vector para colocar todos:
dmop1Ranges = data.frame("ut_ms" = globalMS)

#Completar ventana y graficar:
for (r in 1:nrow(dPairs)) {
  na = !is.na(match(as.character(dPairs[r,1]),colnames(dmop1Scale))) && !is.na(match(as.character(dPairs[r,2]),colnames(dmop1Scale)))
  if(na) {
    dRange = dmop1Scale[, c(1,match(as.character(dPairs[r,1]),colnames(dmop1Scale)),match(as.character(dPairs[r,2]),colnames(dmop1Scale)))]
    name = paste(as.character(dPairs[r,1]), as.character(dPairs[r,2]), sep="-")
    names(dRange)[3] = "end"
    names(dRange)[2] = "start"
    dRange$duration = with(dRange, ave(start-end, 1, FUN = cumsum))
    dRange$duration[dRange$end == 1] = 1
    dRange$duration[dRange$duration == -1] = 1
    names(dRange)[4] = name
    dmop1Ranges = cbind(dmop1Ranges, dRange[4])
  }
}
dmop1RangesMelt <- melt(dmop1Ranges, id='ut_ms')
ggplot(dmop1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

Para estas variables, las ventanas que generan las variaciones son de corto tiempo, por lo que se analizarán posteriormente en una ventana de menor tamaño.

Entrenamiento y Predicción

Promedios por hora, valores para un año:

Para los promedios por hora se tomaron solo los sets de SAAF y LTDATA. Los otros sets se analizarán posteriormente para comparar la correlación que existe entre los circuitos eléctricos y éstos.

Posteriormente se verá que los otros datasets tienen valores que son más cortos en algunos casos que los muestreos de los circuitos de corriente, por lo que se tratarán de forma distinta.

#Cargar los valores
ltdata1 <- read.csv("Data\\train_set\\context--2008-08-22_2010-07-10--ltdata.csv")
saaf1 <- read.csv("Data\\train_set\\context--2008-08-22_2010-07-10--saaf.csv")
power1 <- read.csv("Data\\train_set\\power--2008-08-22_2010-07-10.csv")

#Convertir el tiempo
power1$ut_ms <- as.POSIXct((((power1['ut_ms'])/1000)[,]), origin="1970-01-01")
saaf1$ut_ms <- as.POSIXct((((saaf1['ut_ms'])/1000)[,]), origin="1970-01-01")
ltdata1$ut_ms <- as.POSIXct((((ltdata1['ut_ms'])/1000)[,]), origin="1970-01-01")

Previo al procesamiento de cada set, se agrupan los términos por sus valores por hora:

#Sacar promedios por hora
powerHour = power1
powerHour$ut_ms = cut(powerHour$ut_ms, breaks="hour")
power1HourMean = powerHour %>% group_by(ut_ms) %>% summarise_each(funs(mean))

saaf1Hour = saaf1
saaf1Hour$ut_ms = cut(saaf1Hour$ut_ms, breaks="hour")
saaf1HourMean = saaf1Hour %>% group_by(ut_ms) %>% summarise_each(funs(mean))

ltdata1Hour = ltdata1
ltdata1Hour$ut_ms = cut(ltdata1Hour$ut_ms, breaks="hour")
ltdata1HourMean = ltdata1Hour %>% group_by(ut_ms) %>% summarise_each(funs(mean))

Posteriormente es necesario realizar el Match entre las escalas temporales de los dos sets. Para esto se toma como referencia la escala con más valores, es decir, la de los circuitos eléctricos.

#Match: TimeScale
power1HourMeanMS <- power1HourMean$ut_ms

for (i in 1:nrow(saaf1HourMean)) {
  saaf1HourMean$ut_ms[i] <- power1HourMeanMS[findInterval(saaf1HourMean$ut_ms[i],power1HourMeanMS)]
}

for (i in 1:nrow(ltdata1HourMean)) {
  ltdata1HourMean$ut_ms[i] <- power1HourMeanMS[findInterval(ltdata1HourMean$ut_ms[i],power1HourMeanMS)]
}

Previo al entrenamiento es necesario construir un frame con todas las características que influyen sobre la predicción. En este sentido se tienen los sets de SAAF, y LTDATA.

#Merge:
power1Merge = power1HourMean
power1Merge=merge(x=power1Merge, y=saaf1HourMean, by="ut_ms", all.x=TRUE)
power1Merge=merge(x=power1Merge, y=ltdata1HourMean, by="ut_ms", all.x=TRUE)

Al hacer el merge, varios valores en los datasets agregados tienen valores vacíos. En este caso es necesario completar los valores. Esto se puede realizar mediante interpolación.

power1Merge$sa = na.approx(power1Merge$sa,na.rm = FALSE, rule=2)
power1Merge$sx = na.approx(power1Merge$sx,na.rm = FALSE, rule=2)
power1Merge$sy = na.approx(power1Merge$sy,na.rm = FALSE, rule=2)
power1Merge$sz = na.approx(power1Merge$sz,na.rm = FALSE, rule=2)

power1Merge$sunmars_km = na.spline(power1Merge$sunmars_km, na.rm = FALSE)
power1Merge$earthmars_km = na.spline(power1Merge$earthmars_km, na.rm = FALSE)
power1Merge$sunmarsearthangle_deg = na.spline(power1Merge$sunmarsearthangle_deg, na.rm = FALSE)
power1Merge$solarconstantmars = na.spline(power1Merge$solarconstantmars, na.rm = FALSE)
power1Merge$occultationduration_min = na.spline(power1Merge$occultationduration_min, na.rm = FALSE)
power1Merge$eclipseduration_min = na.approx(power1Merge$eclipseduration_min, na.rm = FALSE, rule=2)

Una vez preparados los datos, se puede proceder a entrenar nuestro modelo.

Random Forest

Se entrenó el RandomForest con nTree = 15. Se entrenan los 33 circuitos con las columnas de los sets SAAF y LTDATA.

#Regresión
rmseSum = 0
i = 1
#for(i in 1:33){
  #Campo a predecir
  predictField = i
  #Columnas en juego:
  predictCols = colnames(power1Merge[,-1])
  #Set de entrenamiento y pruebas
  train = power1Merge[1:12000,-1]
  test = power1Merge[12001:16000,-1]

  colName = predictCols[predictField]

  #Entrenamiento:
  rf = randomForest(as.formula(paste(colName
                                      ,"~ sa + sx + sy + sz + sunmars_km + earthmars_km + sunmarsearthangle_deg + solarconstantmars + eclipseduration_min + occultationduration_min"))
                     ,data=train, ntree=15, proximity=TRUE, importance=TRUE)
  #Prueba:
  predicted = predict(rf, test)
  predCol = test[,c(colName)]
  #Medición error:
  r2 = rmse(predCol, predicted)

  #Graficar Predicción vs. Referencia:
  p = ggplot(aes(x=actual, y=pred), data=data.frame(actual=predCol, pred=predict(rf, test)))
  p = p + geom_point() +  geom_abline(color="red") + ggtitle(paste("RandomForest Regression RMSE=", r2, sep=""))
  #Acumulación Error:
  rmseSum = rmseSum + r2
  #Guardar Imagen
  ggsave(paste("Predict-",i,".png"), p)
#}

Resultados preliminares

errorTotal=rmseSum/33
errorTotal

plot(rf)
print(rf)
importance(rf)
varImpPlot(rf)
MDSplot(rf)

En la gráfica para el modelo entrenado, se puede apreciar que se pueden considerar 15~20 como el número tolerable de árboles para el sistema. Se podría decender más en la cantidad, sin embargo los tiempos de procesamiento serían muy altos.

La predicción realizada tiene valores de error muy bajos. En el mejor de los casos, se pudo predecir con un \(\epsilon\) de \(2.0744e^{-08}\). En el peor de los casos, el error fue de \(0.10185\). Al sumar el error para los 33 circuitos se pudo llegar a un error de \(0.05079785\), lo cual se considera suficiente para esta experiencia.

Caret + gbm

Stochastic Gradient Boosting

Como alternativa a randomforest, se probó con otro package: Caret http://topepo.github.io/caret/training.html#example. Aún utilizando un CrossValidation de 2, no se pudo concluir el entrenamiento del sistema después de 15 minutos de espera. Se deja como referencia para futuras investigaciones.

install.packages("caret")
library(caret)
inTraining <- createDataPartition(power1Merge$NPWD2372, p = .75, list = FALSE)
training <- power1Merge[ inTraining,]
testing  <- power1Merge[-inTraining,]

fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 2,
                           ## repeated ten times
                           repeats = 2)

gbmFit1 <- train(NPWD2372 ~ ., data = training,
                 method = "gbm",
                 trControl = fitControl,
                 verbose = FALSE)

gbmFit1

#}

Caret + svmRadial

Support Vector Machines with Radial Basis Function Kernel
inTraining <- createDataPartition(power1Merge$NPWD2372, p = .75, list = FALSE)
training <- power1Merge[ inTraining,]
testing  <- power1Merge[-inTraining,]

fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 2,
                           ## repeated ten times
                           repeats = 2)

svmFit <- train(NPWD2372 ~ ., data = training,
                 method = "svmRadial",
                 trControl = fitControl,
                 preProcess = c("center", "scale"),
                 tuneLength = 8,
                 metric = "ROC")
svmFit

#}

Caret + rda

Regularized Discriminant Analysis
inTraining <- createDataPartition(power1Merge$NPWD2372, p = .75, list = FALSE)
training <- power1Merge[ inTraining,]
testing  <- power1Merge[-inTraining,]

fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 2,
                           ## repeated ten times
                           repeats = 2)

rdaFit <- train(NPWD2372 ~ ., data = training,
                 method = "rda",
                 trControl = fitControl,
                 tuneLength = 4,
                 metric = "ROC")
rdaFit

#}

Rango menor de tiempo

Una alternativa al proceso realizado anteriormente, considera tomar los tiempos en segundos en la menor escala posible, la entregada por las muestras de potencia. Hacer un match de todos los valores a esta escala.

Debido al alto costo computacional requerido para esto, se ha realizado esto para un rango de un día.

Debido al volumen de los datos, se utilizará una ventana de tiempo, llamada por simplicidad Enero2010. Esta corresponde a un día, desde las 21:00 del 31.12.2009 hasta las 21:00 del 01.01.2010.

power1 = power1[which(power1$ut_ms>"2009-12-31 20:59:00" & power1$ut_ms<"2010-01-01 21:01:00"),]
startDate = power1[1,1]
endDate = power1[nrow(power1),1]

De esta manera, para los demás sets se tiene:

ltdata1 = ltdata1[which(ltdata1$ut_ms>=startDate & ltdata1$ut_ms<=endDate),]
saaf1 = saaf1[which(saaf1$ut_ms>=startDate & saaf1$ut_ms<=endDate),]
dmop1 = dmop1[which(dmop1$ut_ms>=startDate & dmop1$ut_ms<=endDate),]
power1 = power1[which(power1$ut_ms>=startDate & power1$ut_ms<=endDate),]
ftl1 = ftl1[which(ftl1$utb_ms>=startDate & ftl1$ute_ms<=endDate),]
evtf1 = evtf1[which(evtf1$ut_ms>=startDate & evtf1$ut_ms<=endDate),]

Visualización: Enero 01 2010

A nivel de detalle para un día, se pueden apreciar todas las gráficas de corriente para los 33 circuitos en la siguiente gráfica:

power1Melt <- melt(power1, id='ut_ms')
ggplot(power1Melt, aes(x=ut_ms, y=value, color=variable))+geom_line()

En mayor detalle, los 500 primeros valores:

La gráfica de los filtros realizados anteriormente para los distintos sets se puede apreciar para la misma ventana:

globalMS = power1$ut_ms
sum <- rowSums(power1[,c(2:34)])
ut_ms<-power1$ut_ms
power1Sum <- data.frame(ut_ms,sum)
ggplot(power1Sum,aes(x=ut_ms, y=sum)) + geom_line()

Valores SAAF1
#Para los valores originales, la curva de SA:
ggplot(saaf1,aes(x=ut_ms, y=sa)) + geom_line()

#Haciendo el match de los valores:
saaf1Ranges=data.frame("ut_ms" = power1$ut_ms)
saaf1$ut_ms = power1$ut_ms[sapply(saaf1$ut_ms, function(x) findInterval(x, power1$ut_ms, all.inside = TRUE, rightmost.closed = TRUE))]
saaf1Ranges=merge(x=saaf1Ranges, y=saaf1, by="ut_ms", all.x=TRUE)
saaf1Ranges$sa = na.approx(saaf1Ranges$sa,na.rm = FALSE, rule=2)
saaf1Ranges$sx = na.approx(saaf1Ranges$sx,na.rm = FALSE, rule=2)
saaf1Ranges$sy = na.approx(saaf1Ranges$sy,na.rm = FALSE, rule=2)
saaf1Ranges$sz = na.approx(saaf1Ranges$sz,na.rm = FALSE, rule=2)

#Juntar los pares:
saaf1RangesMelt <- melt(saaf1Ranges, id='ut_ms')
ggplot(saaf1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

Valores LTDATA
ltdata1Ranges=data.frame("ut_ms" = power1$ut_ms)
ltdata1$ut_ms = power1$ut_ms[sapply(ltdata1$ut_ms, function(x) findInterval(x, power1$ut_ms, all.inside = TRUE, rightmost.closed = TRUE))]
ltdata1Ranges=merge(x=ltdata1Ranges, y=ltdata1, by="ut_ms", all.x=TRUE)
ltdata1Ranges$sunmars_km = na.spline(ltdata1Ranges$sunmars_km, na.rm = FALSE)
ltdata1Ranges$earthmars_km = na.spline(ltdata1Ranges$earthmars_km, na.rm = FALSE)
ltdata1Ranges$sunmarsearthangle_deg = na.spline(ltdata1Ranges$sunmarsearthangle_deg, na.rm = FALSE)
ltdata1Ranges$solarconstantmars = na.spline(ltdata1Ranges$solarconstantmars, na.rm = FALSE)
ltdata1Ranges$occultationduration_min = na.spline(ltdata1Ranges$occultationduration_min, na.rm = FALSE)
ltdata1Ranges$eclipseduration_min = na.approx(ltdata1Ranges$eclipseduration_min, na.rm = FALSE, rule=2)


#Juntar los pares:
ltdata1RangesMelt <- melt(ltdata1Ranges, id='ut_ms')
ggplot(ltdata1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

Valores EVTF
evtf1Filter = evtf1[grepl("UMBRA", evtf1$description),]
evtf1Scale = data.frame("ut_ms" = globalMS)
evtf1Filter$ut_ms = globalMS[sapply(evtf1Filter$ut_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]
evtf1Groups = dcast(evtf1Filter, ut_ms ~ description, length)
evtf1Scale=merge(x=evtf1Scale, y=evtf1Groups, by="ut_ms", all.x=TRUE)
evtf1Scale[is.na(evtf1Scale)] = 0
ePairs = data.frame(rbind(c("DEI_PENUMBRA_START","DEI_PENUMBRA_END"), 
                          c("MAR_PENUMBRA_START", "MAR_PENUMBRA_END"),
                          c("MAR_UMBRA_START", "MAR_UMBRA_END"),
                          c("PHO_PENUMBRA_START", "PHO_PENUMBRA_END"),
                          c("PHO_UMBRA_START", "PHO_UMBRA_END")))

evtf1Ranges = data.frame("ut_ms" = evtf1Scale$ut_ms)

#Iterar sobre los pares:
for (r in 1:nrow(ePairs)) {
  na = !is.na(match(as.character(ePairs[r,1]),colnames(evtf1Scale))) && !is.na(match(as.character(ePairs[r,2]),colnames(evtf1Scale)))
  if(na){
    eRange = evtf1Scale[, c(1,match(as.character(ePairs[r,1]),colnames(evtf1Scale)),match(as.character(ePairs[r,2]),colnames(evtf1Scale)))]
    names(eRange)[3] = "end"
    names(eRange)[2] = "start"
    eRange$duration = with(eRange, ave(start-end, 1, FUN = cumsum))
    eRange$duration[eRange$end == 1] = 1
    names(eRange)[4] = paste(as.character(ePairs[r,1]), as.character(ePairs[r,2]), sep = " - ")
    evtf1Ranges = cbind(evtf1Ranges, eRange[4])
  }
}

#Juntar los pares:
evtf1RangesMelt <- melt(evtf1Ranges, id='ut_ms')
ggplot(evtf1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

Valores FTL
ftl1Ranges = data.frame("ut_ms" = ftl1Scale$ut_ms)

ftl1Filt = ftl1[with(ftl1, grepl("FALSE", ftl1$flagcomms)&grepl("SLEW", ftl1$type)),]
#Match de los tiempos con el global:
ftl1Filt$utb_ms = globalMS[sapply(ftl1Filt$utb_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]
ftl1Filt$ute_ms = globalMS[sapply(ftl1Filt$ute_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]
#Crear un vector vacío
ftl1Scale = data.frame("ut_ms" = globalMS)
ftl1Scale$type = ""
ftl1Scale$flagcomms = ""
ftl1Scale$start = 0
ftl1Scale$end = 0
#Iterar sobre la colección:
for (i in 1:nrow(ftl1Filt)) {
  x = findInterval(ftl1Filt$utb_ms[i], globalMS, all.inside = TRUE, rightmost.closed = TRUE)
  y = findInterval(ftl1Filt$ute_ms[i], globalMS, all.inside = TRUE, rightmost.closed = TRUE)

  if(ftl1Scale[x,]$start == 1 || ftl1Scale[x,]$end == 1 || ftl1Scale[y,]$start == 1 || ftl1Scale[x,]$end == 1){

  }
  else{
    ftl1Scale[x,]$type = as.character(ftl1Filt$type[i])
    ftl1Scale[x,]$flagcomms = as.character(ftl1Filt$flagcomms[i])
    ftl1Scale[x,]$start = 1

    ftl1Scale[y,]$type = as.character(ftl1Filt$type[i])
    ftl1Scale[y,]$flagcomms = as.character(ftl1Filt$flagcomms[i])
    ftl1Scale[y,]$end = 1
  }
}
fRange = ftl1Scale[, c("ut_ms","start","end")]
fRange$duration = with(fRange, ave(start-end, 1, FUN = cumsum))
fRange$duration[fRange$end == 1] = 1
#ggplot(fRange,aes(x=ut_ms, y=duration)) + geom_line()

ftl1EF = ftl1Scale
fRange = ftl1EF[, c("ut_ms","start","end")]
fRange$duration = with(fRange, ave(start-end, 1, FUN = cumsum))
fRange$duration[fRange$end == 1] = 1
names(fRange)[4] = "EarthFalse"
ftl1Ranges = cbind(ftl1Ranges, fRange[4])

ftl1Filt = ftl1[with(ftl1, grepl("TRUE", ftl1$flagcomms)&grepl("EARTH", ftl1$type)),]
#Match de los tiempos con el global:
ftl1Filt$utb_ms = globalMS[sapply(ftl1Filt$utb_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]
ftl1Filt$ute_ms = globalMS[sapply(ftl1Filt$ute_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]
#Crear un vector vacío
ftl1Scale = data.frame("ut_ms" = globalMS)
ftl1Scale$type = ""
ftl1Scale$flagcomms = ""
ftl1Scale$start = 0
ftl1Scale$end = 0
#Iterar sobre la colección:
for (i in 1:nrow(ftl1Filt)) {
  x = findInterval(ftl1Filt$utb_ms[i], globalMS, all.inside = TRUE, rightmost.closed = TRUE)
  y = findInterval(ftl1Filt$ute_ms[i], globalMS, all.inside = TRUE, rightmost.closed = TRUE)

  if(ftl1Scale[x,]$start == 1 || ftl1Scale[x,]$end == 1 || ftl1Scale[y,]$start == 1 || ftl1Scale[x,]$end == 1){

  }
  else{
    ftl1Scale[x,]$type = as.character(ftl1Filt$type[i])
    ftl1Scale[x,]$flagcomms = as.character(ftl1Filt$flagcomms[i])
    ftl1Scale[x,]$start = 1

    ftl1Scale[y,]$type = as.character(ftl1Filt$type[i])
    ftl1Scale[y,]$flagcomms = as.character(ftl1Filt$flagcomms[i])
    ftl1Scale[y,]$end = 1
  }
}
fRange = ftl1Scale[, c("ut_ms","start","end")]
fRange$duration = with(fRange, ave(start-end, 1, FUN = cumsum))
fRange$duration[fRange$end == 1] = 1

ftl1ET = ftl1Scale
fRange = ftl1ET[, c("ut_ms","start","end")]
fRange$duration = with(fRange, ave(start-end, 1, FUN = cumsum))
fRange$duration[fRange$end == 1] = 1
names(fRange)[4] = "EarthTrue"
ftl1Ranges = cbind(ftl1Ranges, fRange[4])

ftl1RangesMelt <- melt(ftl1Ranges, id='ut_ms')
ggplot(ftl1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

Valores DMOP
dmop1 = dmop1[grepl("^A", dmop1$subsystem),]
dmop1$group = sapply(dmop1$subsystem, substring, 0, 4)
dmop1$end = sapply(dmop1$subsystem, str_sub, -2,-1)
dmop1$status = sapply(dmop1$subsystem, function(x) paste(sprintf("%04d", strtoi(intToBin(str_sub(x,-4,-4)))),sprintf("%04d", strtoi(intToBin(str_sub(x,-3,-3)))),sep = ""))
#Generar un listado de pares:
dPairs = data.frame(rbind(c("ATTTF310A","ATTTF310B"), 
                          c("AAAAF85A1", "AAAAF84A1"),
                          c("AAAAF57A1", "AAAAF60A1"),
                          c("AAAAF59A1", "AAAAF93A1"),
                          c("ASSSF06A0", "ASSSF01A0"),
                          c("AMMMF19A0", "AMMMF18A0"),
                          c("AMMMF24A0", "AMMMF23A0"),
                          c("AAAAF20D1", "AAAAF18D1"),
                          c("AXXX305B", "AXXX301B"),
                          c("AMMMF20A0", "AMMMF11A0")))
names(dPairs)[1] = "ON"
names(dPairs)[2] = "OFF"
#Filtrar aquellos que estén entre los pares:
dmopFilt = dmop1[which(dmop1$subsystem %in% dPairs$ON | dmop1$subsystem %in%  dPairs$OFF),]
dmopFilt$ut_ms = globalMS[sapply(dmopFilt$ut_ms, function(x) findInterval(x, globalMS, all.inside = TRUE, rightmost.closed = TRUE))]
#Transponer los filtrados:
dmopGroups = dcast(dmopFilt, ut_ms ~ subsystem, length)
dmop1Scale = data.frame("ut_ms" = globalMS)
dmop1Scale=merge(x=dmop1Scale, y=dmopGroups, by="ut_ms", all.x=TRUE)
dmop1Scale[is.na(dmop1Scale)] = 0

#Crear vector para colocar todos:
dmop1Ranges = data.frame("ut_ms" = globalMS)
#Completar ventana y graficar:
for (r in 1:nrow(dPairs)) {
  na = !is.na(match(as.character(dPairs[r,1]),colnames(dmop1Scale))) && !is.na(match(as.character(dPairs[r,2]),colnames(dmop1Scale)))
  if(na) {
    dRange = dmop1Scale[, c(1,match(as.character(dPairs[r,1]),colnames(dmop1Scale)),match(as.character(dPairs[r,2]),colnames(dmop1Scale)))]
    name = paste(as.character(dPairs[r,1]), as.character(dPairs[r,2]), sep="-")
    names(dRange)[3] = "end"
    names(dRange)[2] = "start"
    dRange$duration = with(dRange, ave(start-end, 1, FUN = cumsum))
    dRange$duration[dRange$end == 1] = 1
    dRange$duration[dRange$duration == -1] = 1
    names(dRange)[4] = name
    dmop1Ranges = cbind(dmop1Ranges, dRange[4])
  }
}
dmop1RangesMelt <- melt(dmop1Ranges, id='ut_ms')
ggplot(dmop1RangesMelt, aes(x=ut_ms, y=value, color=variable))+geom_line()

Correlaciones

Para estos valores no se realizará predicción, ya que solo se consideran unos pocos valores para cada variable. De todas formas, se busca encontrar si existe alguna alta correlación entre los distintos valores.

Para el caso de la ventana, se procedió a unir todos los vectores y buscar construir una matriz de correlaciones. Al juntar todos los vectores, la matriz tiene 55 variables sin contar la escala temporal. La construcción de una matriz de esta forma es relativamente complicado, por lo que procederá a trabajar con uniones más pequeñas.

power1Merge = power1
power1Merge=merge(x=power1Merge, y=dmop1Ranges, by="ut_ms", all.x=TRUE)
power1Merge=merge(x=power1Merge, y=evtf1Ranges, by="ut_ms", all.x=TRUE)
power1Merge=merge(x=power1Merge, y=ftl1Ranges, by="ut_ms", all.x=TRUE)
power1Merge=merge(x=power1Merge, y=saaf1Ranges, by="ut_ms", all.x=TRUE)
power1Merge=merge(x=power1Merge, y=ltdata1Ranges, by="ut_ms", all.x=TRUE)
power1Merge[is.na(power1Merge)] <- 0
write.csv(power1Merge, "Web/Data/Power1MergeWindow.csv", row.names = FALSE)

El resultado de obtener una matriz de correlaciones con tantos valores genera mucho ruido. Se procederá a calcular en pares las relaciones contra los circuitos de poder, omitiendo que puedan existir entre las otras variables.

POWER y DMOP

power1Merge = power1
power1Merge=merge(x=power1Merge, y=dmop1Ranges, by="ut_ms", all.x=TRUE)
power1Merge[is.na(power1Merge)] <- 0
#Es necesario eliminar la columna de tiempo
power1Merge[1] = NULL
library(corrplot)
corrplot(cor(power1Merge), method = "square", tl.pos="n")

Se puede apreciar que existe una baja correlación entre los comandos elegidos en DMOP y los valores de corriente de los circuitos. Se podría intentar con otros pares de comandos.

POWER y EVTF

power1Merge = power1
power1Merge=merge(x=power1Merge, y=evtf1Ranges, by="ut_ms", all.x=TRUE)
power1Merge[is.na(power1Merge)] <- 0
#Es necesario eliminar la columna de tiempo
power1Merge[1] = NULL
corrplot(cor(power1Merge), method = "circle", tl.pos="n")

POWER y FTL

power1Merge = power1
power1Merge=merge(x=power1Merge, y=ftl1Ranges, by="ut_ms", all.x=TRUE)
power1Merge[is.na(power1Merge)] <- 0
#Es necesario eliminar la columna de tiempo
power1Merge[1] = NULL

corrplot(cor(power1Merge), method = "square", tl.pos="n")

POWER y SAAF

power1Merge = power1
power1Merge=merge(x=power1Merge, y=saaf1Ranges, by="ut_ms", all.x=TRUE)
power1Merge[is.na(power1Merge)] <- 0
#Es necesario eliminar la columna de tiempo
power1Merge[1] = NULL

corrplot(cor(power1Merge), method = "square", tl.pos="n")

POWER y LTDATA

power1Merge = power1
power1Merge=merge(x=power1Merge, y=ltdata1Ranges, by="ut_ms", all.x=TRUE)
power1Merge[is.na(power1Merge)] <- 0
#Es necesario eliminar la columna de tiempo
power1Merge[1] = NULL

corrplot(cor(power1Merge), method = "square", tl.pos="n")

Conclusiones

Durante el (pre)procesamiento de los datos se encontraron varias dificultades, principalmente relacionadas con las capacidades computaciones. Procesar todos los datos sin agruparlos por hora conllevó tiempos extendidos, este tradeoff seguramente podría incrementar los valores de predicción, sin embargo se descarta para el desarrollo de este trabajo.

De los cálculos realizados se desprende como primera observación que el consumo de potencia tiene una estrecha relación con la cantidad de energía generada, lo que es intuitivo. Esta idea afirma la hipótesis de que existe una predicción gruesa y una fina.

Para poder considerar los otros sets de datos se considera realizar un pre-procesamiento adicional a los datos, como por ejemplo, encontrar la relación entre el consumo de los circuitos y la emisión de los comandos. De la misma manera, la incidencia solar tiene una dependencia sobre los eclipses solares.

Si bien se probaron varios pares de comandos, no se pudo encontrar una correlación entre éstos y los circuitos. Tampoco se pudo probar con extensión estos ya que los cálculos se realizaron para las mediciones de un día particular. Lo que no puede ser tomado como referencia.

Las correlaciones encontradas dependen de los mismos resultados anteriores, por lo que no se pueden tomar como referencia.

Pese a los resultados anteriores, se lograron buenos resultados sobre la predicción del consumo eléctrico utilizando RandomForest de 15 árboles. El RMSE se contabilizó cercano al 6%, lo que es superior a lo que se podría esperar para una primera experiencia.

Lo anterior también corresponde a que los set de datos provenían de una fuente límpia. Si bien no tenían una descripción sólida sobre la funcionalidad de los distintos comandos y mensajes, estaban completos, sin vacíos y sin datos incoherentes o outliers.

Finalmente se probó con otras librerías (distintas a randomforest). Se probó con Caret: svmRadial (Support Vector Machines with Radial Basis Function Kernel),rda (Regularized Discriminant Analysis) y gbm (Stochastic Gradient Boosting). En ninguno de los otros casos se pudo llegar a resultados. Se considera que se debe principalmente al volumen de datos que se debían procesar.

Si bien existen otras alternativas para el procesamiento de datos, como ff, big.tables o servicios en la nube, no se probó con estas alternativas ya que el grueso del trabajo realizado consistió en ordenar los datos, encontrar correlaciones, entrenar el modelo y aprender la metodología de minería de datos.