This R Notebook utilizes the javascript functions from the CBE Comfort Tool, freely available on the website’s Github in order to estimate the thermal comfort of the students of KWest, KWEast and Illima projects.
The data used here is associated to the ERDL frog3 database, which for the purposes of analysis exists in the sensor_info2 and all_sensor_data tables.
s <- suppressPackageStartupMessages
library(data.table)
library(lubridate)
library(knitr)
sensors <- fread("data/sensor_info2.csv")
readings <- fread("data/all_sensor_data.csv",showProgress=FALSE)The sensor_info2 table contains various information associated to the reading:
kable(sensors)| sensor_id | building | location | type | end_use |
|---|---|---|---|---|
| 1 | KWWest | Ceiling | Illuminance_ft-c | |
| 2 | KWWest | W wall | Illuminance_ft-c | |
| 3 | KWWest | room | RH_pct | |
| 4 | KWWest | W wall | Temp_F_air | |
| 5 | KWWest | E wall | Temp_F_air | |
| 6 | KWWest | W wall | Temp_F_surface | |
| 7 | KWWest | E wall | Temp_F_surface | |
| 8 | KWEast | Ceiling | Illuminance_ft-c | |
| 9 | KWEast | W wall | Illuminance_ft-c | |
| 10 | KWEast | room | RH_pct | |
| 11 | KWEast | W wall | Temp_F_air | |
| 12 | KWEast | E wall | Temp_F_air | |
| 13 | KWEast | W wall | Temp_F_surface | |
| 14 | KWEast | E wall | Temp_F_surface | |
| 15 | Illima | Ceiling | Illuminance_ft-c | |
| 16 | Illima | W Wall | Illuminance_ft-c | |
| 17 | Illima | E wall | Temp_F_air | |
| 18 | Illima | W wall | Temp_F_air | |
| 19 | Illima | W wall | Temp_F_surface | |
| 20 | Illima | E wall | Temp_F_surface | |
| 21 | Illima | room | RH_pct | |
| 22 | EwaWeather | RH_pct | ||
| 23 | EwaWeather | Solar_Radiation_Wm-2 | ||
| 24 | EwaWeather | Temp_F_air | ||
| 25 | KawaikiniWeather | RH_pct | ||
| 26 | KawaikiniWeather | Solar_Radiation_Wm-2 | ||
| 27 | KawaikiniWeather | Temp_F_air | ||
| 28 | Illima | Energy_kwh | ceiling_fans | |
| 29 | Illima | Energy_kwh | condensing_unit | |
| 30 | Illima | Energy_kwh | exhaust_fans | |
| 31 | Illima | Energy_kwh | fan_coil_unit | |
| 32 | Illima | Energy_kwh | lighting_exterior | |
| 33 | Illima | Energy_kwh | lighting_main_space | |
| 34 | Illima | Energy_kwh | panel_feed | |
| 35 | Illima | Energy_kwh | pv_inverter_1 | |
| 36 | Illima | Energy_kwh | pv_inverter_2 | |
| 41 | KWEast | Energy_kwh | ceiling_fans | |
| 42 | KWEast | Energy_kwh | condensing_unit | |
| 43 | KWEast | Energy_kwh | exhaust_fans | |
| 44 | KWEast | Energy_kwh | fan_coil_unit | |
| 45 | KWEast | Energy_kwh | lighting_exterior | |
| 46 | KWEast | Energy_kwh | lighting_main_space | |
| 47 | KWEast | Energy_kwh | panel_feed | |
| 48 | KWEast | Energy_kwh | plug_load | |
| 49 | KWEast | Energy_kwh | louver_actuator | |
| 50 | KWEast | Energy_kwh | pv_generation | |
| 51 | KWWest | Energy_kwh | ceiling_fans | |
| 52 | KWWest | Energy_kwh | condensing_unit | |
| 53 | KWWest | Energy_kwh | exhaust_fans | |
| 54 | KWWest | Energy_kwh | fan_coil_unit | |
| 55 | KWWest | Energy_kwh | lighting_exterior | |
| 56 | KWWest | Energy_kwh | lighting_main_space | |
| 57 | KWWest | Energy_kwh | panel_feed | |
| 58 | KWWest | Energy_kwh | plug_load | |
| 59 | KWWest | Energy_kwh | louver_actuator | |
| 60 | KWWest | Energy_kwh | pv_generation | |
| 37 | Illima | Energy_kwh | plug_load | |
| 38 | Illima | Energy_kwh | louver_actuator | |
| 39 | Illima | Energy_kwh | pv_generation |
While all_sensor_data contain the readings, as shown on the following sample:
readings <- readings[,.(date=ymd_hms(datetime),sensor_id,reading)]
kable(head(readings))| date | sensor_id | reading |
|---|---|---|
| 2013-08-17 22:15:00 | 49 | 0.00 |
| 2013-08-17 22:25:00 | 41 | 0.00 |
| 2013-08-17 22:25:00 | 42 | 0.00 |
| 2013-08-17 22:25:00 | 43 | 0.00 |
| 2013-08-17 22:25:00 | 44 | 0.00 |
| 2013-08-17 22:25:00 | 45 | 0.04 |
To calculate the thermal comfort, we first separate the readings associated to each building.
dt <- merge(sensors,readings,by="sensor_id")
west <- dt[building=="KWWest"]
east <- dt[building=="KWEast"]
illima <- dt[building=="Illima"]Next, we format each building table readings, so that each sensor is in a specific column. Among the required sensor readings for the thermal comfort model, we currently only have absolute temperature and relative humidity, and rows which have missing either temperature or relative humidity are filtered as the thermal comfort model does not handle incomplete observations:
LongFormatReadings <- function(building){
building <- dcast(building,date ~ sensor_id,value.var=c("reading"))
building <- building[,c("date","4","3"),with=FALSE]
building <- building[complete.cases(building)]
colnames(building) <- c("date","temp","relh")
return(building)
}
west <- LongFormatReadings(west)
#east <- LongFormatReadings(east)
#illima <- LongFormatReadings(illima)For the remaining required variables, we made the following assumptions given the setting of the 3 buildings were a classroom:
AddOtherVariableAssumptions <- function(building){
building$radtemp <- building$temp
building$vel <- 0.1
building$met <- 1.1
building$clo <- 0.5
building$wme <- 0
return(building)
}
west <- AddOtherVariableAssumptions(west)As shown below (sample):
kable(head(west))| date | temp | relh | radtemp | vel | met | clo | wme |
|---|---|---|---|---|---|---|---|
| 2014-07-01 00:05:00 | 81.97244 | 57.73286 | 81.97244 | 0.1 | 1.1 | 0.5 | 0 |
| 2014-07-01 00:15:00 | 81.84932 | 57.11775 | 81.84932 | 0.1 | 1.1 | 0.5 | 0 |
| 2014-07-01 00:25:00 | 81.84932 | 59.05097 | 81.84932 | 0.1 | 1.1 | 0.5 | 0 |
| 2014-07-01 00:35:00 | 81.72619 | 59.35852 | 81.72619 | 0.1 | 1.1 | 0.5 | 0 |
| 2014-07-01 00:45:00 | 81.72619 | 59.35852 | 81.72619 | 0.1 | 1.1 | 0.5 | 0 |
| 2014-07-01 00:55:00 | 81.60307 | 59.35852 | 81.60307 | 0.1 | 1.1 | 0.5 | 0 |
We use the functions from CBE Comfort Tool website on the pre-processed data to obtain the thermal comfort. The R Library (V8) is used to execute javascript, and return a dataframe object. This Notebook uses 2 functions:
The PMV Elevated Air function is a superset of the PMV function. When airspeed is beyond 0.15 m/s, then the original PMV model can no longer be used. In these cases, the PMV Elevated Air function performs some pre-processing before calling the PMV model. The intuition behind it is that the SET index is used on the data with elevated air speed, and adjust the sensor readings through temperature and relative humidity to reflect the same thermal comfort situation on a lower air speed, thus comforming to the original PMV model [1].
The Elevated Air speed superset function was proposed by the CBE Comfort Tool group,
#Load Google's Javscript Engine V8 (See https://cran.r-project.org/web/packages/V8/vignettes/v8_intro.html)
library(V8)
#Create a new context
ct <- v8()
#Load Javascript Library for forEach function
ct$source(system.file("js/underscore.js", package="V8"))
#Load local comfortModel javscript library (only modified the path of the libraries)
ct$source("comfortmodels.js")## [1] "function (ta, tr, runningMean, vel) {\n var to = (ta + tr) / 2;\n var coolingEffect = 0;\n if (vel >= 0.2 && to > 25) {\n // calculate cooling effect of elevated air speed\n // when top > 25 degC.\n var coolingEffect = 1.7856 * Math.log(vel) + 2.9835;\n }\n var tComf = 0.33 * runningMean + 18.8;\n if(runningMean > 15){\n var tComfILower = tComf - 2;\n var tComfIUpper = tComf + 2 + coolingEffect;\n var tComfIILower = tComf - 3;\n var tComfIIUpper = tComf + 3 + coolingEffect;\n var tComfIIILower = tComf - 4;\n var tComfIIIUpper = tComf + 4 + coolingEffect;\n } else if (12.73 < runningMean && runningMean < 15){\n var tComfLow = 0.33 * 15 + 18.8;\n var tComfILower = tComfLow - 2;\n var tComfIUpper = tComf + 2 + coolingEffect;\n var tComfIILower = tComfLow - 3;\n var tComfIIUpper = tComf + 3 + coolingEffect;\n var tComfIIILower = tComfLow - 4;\n var tComfIIIUpper = tComf + 4 + coolingEffect;\n } else {\n var tComfLow = 0.33 * 15 + 18.8;\n var tComfILower = tComfLow - 2;\n var tComfIUpper = tComf + 2;\n var tComfIILower = tComfLow - 3;\n var tComfIIUpper = tComf + 3 + coolingEffect;\n var tComfIIILower = tComfLow - 4;\n var tComfIIIUpper = tComf + 4 + coolingEffect;\n }\n var acceptabilityI, acceptabilityII, acceptabilityIII;\n\n if (comf.between(to, tComfILower, tComfIUpper)) {\n // compliance at all levels\n acceptabilityI = acceptabilityII = acceptabilityIII = true;\n } else if (comf.between(to, tComfIILower, tComfIIUpper)) {\n // compliance at II and III only\n acceptabilityII = acceptabilityIII = true;\n acceptabilityI = false;\n } else if (comf.between(to, tComfIIILower, tComfIIIUpper)) {\n // compliance at III only\n acceptabilityIII = true;\n acceptabilityI = acceptabilityII = false;\n } else {\n // neither\n acceptabilityI = acceptabilityII = acceptabilityIII = false;\n }\n r = {}\n r.acceptabilityI = acceptabilityI\n r.acceptabilityII = acceptabilityII\n r.acceptabilityIII = acceptabilityIII\n r.tComfILower = tComfILower\n r.tComfIILower = tComfIILower\n r.tComfIIILower = tComfIIILower\n r.tComfIUpper = tComfIUpper\n r.tComfIIUpper = tComfIIUpper\n r.tComfIIIUpper = tComfIIIUpper\n return r\n return [[acceptabilityIII, tComfIIILower, tComfIIIUpper],\n [acceptabilityII, tComfIILower, tComfIIUpper],\n [acceptabilityI, tComfILower, tComfIUpper]];\n}"ct$source("util.js")## [1] "function (x) {\n return (x - 32) * 5 / 9;\n}"ct$source("psychrometrics.js")## [1] "function (ta, vel, tglobe, diameter, emissivity) {\n pow = Math.pow;\n return pow(pow(tglobe + 273, 4) + (1.1 * pow(10, 8) * pow(vel, 0.6)) / (emissivity * pow(diameter, 0.4)) * (tglobe - ta), 0.25) - 273;\n}"#Apply the function over all the table for pmvElevatedAirspeed
west.elevated.air.pmv <- data.table(ct$call("_.map", west, JS("function(x){return(comf.pmvElevatedAirspeed(util.FtoC(x.temp),util.FtoC(x.radtemp),x.vel/196.85,x.relh,x.clo,x.met,x.wme))}")))
west.pmv<- data.table(ct$call("_.map", west, JS("function(x){return(comf.pmv(util.FtoC(x.temp),util.FtoC(x.radtemp),x.vel/196.85,x.relh,x.clo,x.met,x.wme))}")))The CBE Comfort Tool functions only returns the PMV scores. Depending on the Standard adopted, the threshold of the PMV score varies in indicating wether the group of people of the readings, subject to the PPD value is comfortable or not on the thermal conditions. For ASHRAE55, this threshold occurs on -0.5 and +0.5 (Refer to the standard for a finer granularity comfort label such as Hot, Cold, Very Hot, Very Cold, etc.)
west.elevated.air.pmv$iscomfortable <- ifelse(west.elevated.air.pmv$pmv > -0.5 & west.elevated.air.pmv$pmv < 0.5,"comfortable","uncomfortable")
kable(west.elevated.air.pmv[,.(frequency=length(pmv)),by="iscomfortable"])| iscomfortable | frequency |
|---|---|
| uncomfortable | 53345 |
| comfortable | 763 |
library(dygraphs)
west.elevated.air.pmv$date <- west$date
dygraph(west.elevated.air.pmv[,.(date,pmv)], main = "Comfort Level") %>% dyRangeSelector()[1] Web application for thermal comfort visualization and calculation according to ASHRAE Standard 55.