commit b0f5903d3345be34cb8288e6c9a06945f057931e
parent 37f99067dbc862580be7ca3355ad81e90b9fb283
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Thu, 1 Apr 2021 16:30:22 +0300
Velocity distribution graph.
Diffstat:
13 files changed, 498 insertions(+), 144 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -1,8 +1,8 @@
/build
Rplots.pdf
-# Created by https://www.toptal.com/developers/gitignore/api/vim,linux
-# Edit at https://www.toptal.com/developers/gitignore?templates=vim,linux
+# Created by https://www.toptal.com/developers/gitignore/api/vim,linux,r
+# Edit at https://www.toptal.com/developers/gitignore?templates=vim,linux,r
### Linux ###
*~
@@ -19,6 +19,51 @@ Rplots.pdf
# .nfs files are created when an open file is removed but is still being accessed
.nfs*
+### R ###
+# History files
+.Rhistory
+.Rapp.history
+
+# Session Data files
+.RData
+
+# User-specific files
+.Ruserdata
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user/
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# R Environment Variables
+.Renviron
+
+### R.Bookdown Stack ###
+# R package: bookdown caching files
+/*_files/
+
### Vim ###
# Swap
[._]*.s[a-v][a-z]
@@ -39,4 +84,4 @@ tags
# Persistent undo
[._]*.un~
-# End of https://www.toptal.com/developers/gitignore/api/vim,linux
+# End of https://www.toptal.com/developers/gitignore/api/vim,linux,r
diff --git a/Makefile b/Makefile
@@ -16,6 +16,8 @@ all: build/iccsa-21-wind-slides.pdf
build/main.pdf: main.tex
build/main.pdf: build/inkscape/anemometer.eps
build/main.pdf: build/daily-stats.eps
+build/main.pdf: build/daily-rmse.eps
+build/main.pdf: build/gnuplot/velocity-dist.eps
#build/main.pdf: build/gnuplot/anemometer.eps
#build/main.pdf: build/gnuplot/anemometer.svg
build/main.pdf:
@@ -46,6 +48,9 @@ build/inkscape/%.eps: inkscape/%.svg
build/daily-stats.svg: gnuplot/daily-stats.gnuplot build/daily-stats
gnuplot -d $<
+build/daily-rmse.svg: gnuplot/daily-rmse.gnuplot build/daily-rmse
+ gnuplot -d $<
+
#build/main.zip: build/gnuplot/*.eps main.tex
build/main.zip: main.tex
@mkdir -p build
diff --git a/R/anal.R b/R/anal.R
@@ -2,11 +2,13 @@ library(DBI)
library(fitdistrplus)
library(plotrix)
library(sn)
+library(doParallel)
+library(foreach)
+registerDoParallel(cores=detectCores())
-datetime_format <- "%FT%H:%M"
+source("R/common.R")
-# https://stackoverflow.com/questions/13023274/how-to-do-printf-in-r
-printf <- function(...) cat(sprintf(...))
+datetime_format <- "%FT%H:%M"
# https://stackoverflow.com/questions/32370485/convert-radians-to-degree-degree-to-radians
rad2deg <- function(rad) {(rad * 180) / (pi)}
@@ -31,6 +33,10 @@ power <- function(x,y) {
#abs(x)**y
}
+rmse <- function(estimated, actual){
+ sqrt(mean((estimated - actual)**2))
+}
+
dweibull2 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g=0) {
#ifelse(x<g,
# a1* abs(c1*b1)*(abs(b1*(x-g))**(c1-1))* exp(-(abs(b1*(x-g))**c1)),
@@ -56,9 +62,6 @@ dgumbel <- function (x, location=0, scale=1) {
(1/scale) * exp(-(z + exp(-z)))
}
-drayleigh <- function (x, location=0, scale=1) {
-}
-
# https://www.tutorialspoint.com/r/r_mean_median_mode.htm
# https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode
getmode <- function(v) {
@@ -209,6 +212,7 @@ fit_velocity_pdf <- function (xx, x, density, sign=-1, plot=TRUE, axis="x", dist
})
pdf_x_est <- c[[2]]*dweibull(sign*pdf_x$x,scale=c[[1]],shape=2)
} else {
+ x_mode <- find_max_x(pdf_x$x,pdf_x$density)
c <- tryCatch(
{
model <- nls(pdf_x$density ~ c*dweibull(sign*pdf_x$x,scale=a,shape=b),
@@ -216,7 +220,8 @@ fit_velocity_pdf <- function (xx, x, density, sign=-1, plot=TRUE, axis="x", dist
start=list(a=1,b=1,c=1),
control=list(warnOnly=TRUE),
algorithm="port",
- lower=c(0.1,0.1,0), upper=c(1000,1000,1000))
+ lower=c(0.1,0.1,0),
+ upper=c(1000,1000,1000))
c = summary(model)$coefficients
},
error=function(cond) {
@@ -239,19 +244,20 @@ fit_velocity_pdf <- function (xx, x, density, sign=-1, plot=TRUE, axis="x", dist
list(coefficients=c,
actual_density=pdf_x$density,
estimated_density=pdf_x_est,
- x=pdf_x$x)
+ x=pdf_x$x,
+ residual=max(abs(pdf_x$density-pdf_x_est)))
}
}
-fit_velocity_pdf_bilateral <- function (velocity, axis="x") {
+fit_velocity_pdf_bilateral <- function (velocity, axis="x", timestamp=0, dist="weibull") {
velocity_hist <- normalize_hist(hist(velocity, plot=FALSE, breaks=100))
x = velocity_hist$x;
v = velocity_hist$density;
d = density(velocity,adjust=0.5)
dist_left <- fit_velocity_pdf(x[x<=0], v[x<=0], density=d$y[d$x<=0],
- sign=-1, plot=FALSE, axis=axis, dist="weibull")
+ sign=-1, plot=FALSE, axis=axis, dist=dist)
dist_right <- fit_velocity_pdf(x[x>=0], v[x>=0], density=d$y[d$x>=0],
- sign=1, plot=FALSE, axis=axis, dist="weibull")
+ sign=1, plot=FALSE, axis=axis, dist=dist)
#x <- c(dist_left$x, dist_right$x)
actual_density <- c(dist_left$actual_density, dist_right$actual_density)
estimated_density <- c(dist_left$estimated_density, dist_right$estimated_density)
@@ -261,7 +267,17 @@ fit_velocity_pdf_bilateral <- function (velocity, axis="x") {
#lines(x, estimated_density, col='red', lwd=2)
lines(dist_left$x, dist_left$estimated_density, col='red', lwd=2)
lines(dist_right$x, dist_right$estimated_density, col='blue', lwd=2)
- dist_left$coefficients
+ #printf("RMSE %f\n", rmse(estimated_density, actual_density))
+ path <- file.path('build', 'velocity', axis, dist)
+ dir.create(path, recursive=TRUE)
+ filename <- file.path(path, sprintf("%d", timestamp))
+ write.table(data.frame(x=x,actual=actual_density,estimated=estimated_density),
+ filename, row.names=FALSE, quote=FALSE)
+ write("\n",file=filename,append=TRUE)
+ write.table(data.frame(x=d$x,density=d$y), filename,
+ row.names=FALSE, quote=FALSE, append=TRUE)
+ list(coefficients=dist_left$coefficients,
+ rmse=rmse(estimated_density, actual_density))
}
fit_velocity_pdf_v2 <- function (x, sign=-1, plot=TRUE, axis="x") {
@@ -316,50 +332,55 @@ fit_direction <- function (x, x_mode, plot=TRUE) {
c
}
-mydb <- dbConnect(RSQLite::SQLite(), "samples/load-cell.sqlite3")
-
time_start = time_string_to_timestamp("2021-03-06T00:00")
time_end = time_string_to_timestamp("2021-03-06T23:00")
+#time_start = time_string_to_timestamp("2021-02-21T00:00")
+#time_end = time_string_to_timestamp("2021-03-31T23:00")
time_delta = 60*60*2
print(time_start)
print(time_end)
-cumulative_data <- data.frame(speed=numeric(0),direction=numeric(0),
- x_mode=numeric(0),abs_x_mode=numeric(0),x_sd=numeric(0),
- y_mean=numeric(0),
- acf_x_a=numeric(0),acf_x_b=numeric(0),acf_x_c=numeric(0),
- x_a=numeric(0),x_b=numeric(0),x_c=numeric(0),
- x_d=numeric(0),x_e=numeric(0),x_f=numeric(0),
- # Y
- y_mode=numeric(0),
- y_a=numeric(0),y_b=numeric(0),y_c=numeric(0),
- y_d=numeric(0),y_e=numeric(0),y_f=numeric(0),
- acf_y_a=numeric(0),acf_y_b=numeric(0),acf_y_c=numeric(0))
timestamps <- seq(time_start, time_end, time_delta)
-pb <- txtProgressBar(min = 0, max = length(timestamps), style = 3)
-
-i <- 0
-for (timestamp in timestamps) {
- velocity = dbGetQuery(mydb,
+#pb <- txtProgressBar(min = 0, max = length(timestamps), style = 3)
+
+#i <- 0
+#for (timestamp in timestamps) {
+cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='remove') %dopar% {
+ db <- dbConnect(RSQLite::SQLite(), "samples/load-cell.sqlite3")
+ cumulative_data <- data.frame(speed=numeric(0),direction=numeric(0),
+ x_mode=numeric(0),abs_x_mode=numeric(0),x_sd=numeric(0),
+ y_mean=numeric(0),
+ acf_x_a=numeric(0),acf_x_b=numeric(0),acf_x_c=numeric(0),
+ x_a=numeric(0),x_b=numeric(0),x_c=numeric(0),
+ x_d=numeric(0),x_e=numeric(0),x_f=numeric(0),
+ # Y
+ y_mode=numeric(0),
+ y_a=numeric(0),y_b=numeric(0),y_c=numeric(0),
+ y_d=numeric(0),y_e=numeric(0),y_f=numeric(0),
+ acf_y_a=numeric(0),acf_y_b=numeric(0),acf_y_c=numeric(0),
+ timestamp=numeric(0),
+ x_rmse=numeric(0), y_rmse=numeric(0), z_rmse=numeric(0))
+ velocity = dbGetQuery(db,
"SELECT x,y,z FROM samples WHERE timestamp BETWEEN :t0 AND (:t0 + :dt) ORDER BY timestamp",
params = list(t0 = timestamp, dt=time_delta))
- velocity = dbGetQuery(mydb,
+ velocity = dbGetQuery(db,
"SELECT x,y,z FROM samples
WHERE timestamp BETWEEN :t0 AND (:t0 + :dt)
- AND ABS(x-:xmean)<10000 AND ABS(y-:ymean)<10000 AND ABS(z-:zmean)<10000
+ AND ABS(x-:xmean)<10000 AND ABS(y-:ymean)<10000 AND ABS(z-:zmean)<10000
ORDER BY timestamp",
params = list(t0 = timestamp, dt=time_delta,
xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z)))
+ if (nrow(velocity) == 0) { next }
+ filename <- sprintf("build/rplot-%d.pdf", timestamp)
+ print(filename)
+ pdf(filename)
old_mean_x <- max(abs(velocity$x)) - mean(velocity$x)
- velocity$x = (velocity$x - mean(velocity$x))
- velocity$y = (velocity$y - mean(velocity$y))
- velocity$z = (velocity$z - mean(velocity$z))
- # convert energy to velocity (U^2 -> U)
- velocity$x = sign(velocity$x)*sqrt(abs(velocity$x))
- velocity$y = sign(velocity$y)*sqrt(abs(velocity$y))
- velocity$z = sign(velocity$z)*sqrt(abs(velocity$z))
+ coefficients <- calibration_coefficients()
+ velocity$x = sample_to_speed(velocity$x, coefficients$x1, coefficients$x2)
+ velocity$y = sample_to_speed(velocity$y, coefficients$y1, coefficients$y2)
+ velocity$z = sample_to_speed(velocity$z, coefficients$z1, coefficients$z2)
x_mode <- getmode(velocity$x)
y_mode <- getmode(velocity$y)
#direction_mode <- getmode(velocity$direction)
@@ -367,6 +388,7 @@ for (timestamp in timestamps) {
velocity$speed = speed(velocity$x,velocity$y,velocity$z)
velocity$direction = direction(velocity$x,velocity$y)
new_row = nrow(cumulative_data)+1
+ cumulative_data[new_row, "timestamp"] = timestamp
cumulative_data[new_row, "x_mode"] = x_mode
cumulative_data[new_row, "y_mode"] = y_mode
cumulative_data[new_row, "speed"] = mean(velocity$speed)
@@ -378,23 +400,31 @@ for (timestamp in timestamps) {
#fit <- fitdist(velocity$x, "gumbel", start=list(a=1, b=1))
#plot(fit)
# Velocity X
- x_coefs <- fit_velocity_pdf_bilateral(velocity$x, axis="x")
- cumulative_data[new_row, "x_a"] = x_coefs[[1]]
- cumulative_data[new_row, "x_b"] = x_coefs[[2]]
+ #for (dist in c("weibull2")) {
+ # fit_velocity_pdf_bilateral(velocity$x, axis="x", timestamp=timestamp, dist="weibull2")
+ # fit_velocity_pdf_bilateral(velocity$y, axis="y", timestamp=timestamp, dist="weibull2")
+ # fit_velocity_pdf_bilateral(velocity$z, axis="z", timestamp=timestamp, dist="weibull2")
+ #}
+ x_model <- fit_velocity_pdf_bilateral(velocity$x, axis="x", timestamp=timestamp, dist="weibull")
+ cumulative_data[new_row, "x_a"] = x_model$coefficients[[1]]
+ cumulative_data[new_row, "x_b"] = x_model$coefficients[[2]]
+ cumulative_data[new_row, "x_rmse"] = x_model$rmse
#cumulative_data[new_row, "x_c"] = x_coefs[[3]]
#cumulative_data[new_row, "x_d"] = x_coefs[[4]]
#cumulative_data[new_row, "x_e"] = x_coefs[[5]]
#cumulative_data[new_row, "x_f"] = x_coefs[[6]]
# Velocity Y
- y_coefs <- fit_velocity_pdf_bilateral(velocity$y, axis="y")
- cumulative_data[new_row, "y_a"] = y_coefs[[1]]
- cumulative_data[new_row, "y_b"] = y_coefs[[2]]
+ y_model <- fit_velocity_pdf_bilateral(velocity$y, axis="y", timestamp=timestamp, dist="weibull")
+ cumulative_data[new_row, "y_a"] = y_model$coefficients[[1]]
+ cumulative_data[new_row, "y_b"] = y_model$coefficients[[2]]
+ cumulative_data[new_row, "y_rmse"] = y_model$rmse
#cumulative_data[new_row, "y_c"] = y_coefs[[3]]
#cumulative_data[new_row, "y_d"] = y_coefs[[4]]
#cumulative_data[new_row, "y_e"] = y_coefs[[5]]
#cumulative_data[new_row, "y_f"] = y_coefs[[6]]
# Velocity Z
- z_coefs <- fit_velocity_pdf_bilateral(velocity$z, axis="z")
+ z_model <- fit_velocity_pdf_bilateral(velocity$z, axis="z", timestamp=timestamp, dist="weibull")
+ cumulative_data[new_row, "z_rmse"] = z_model$rmse
# ACF X
x_coefs <- fit_velocity_acf(velocity$x, plot=TRUE, axis="x")
cumulative_data[new_row, "acf_x_a"] = x_coefs[[1]]
@@ -406,7 +436,7 @@ for (timestamp in timestamps) {
cumulative_data[new_row, "acf_y_c"] = y_coefs[[3]]
direction_coefs <- fit_direction(velocity$direction, direction_mode, plot=TRUE)
- print(direction_coefs)
+ #print(direction_coefs)
old_mar = par("mar")
direction_hist <- hist(velocity$direction, plot=FALSE, breaks=100)
@@ -418,19 +448,24 @@ for (timestamp in timestamps) {
cumulative_data[new_row, "abs_x_mode"] = abs(x_mode)
cumulative_data[new_row, "x_sd"] = sd(velocity$x)
cumulative_data[new_row, "y_mean"] = mean(velocity$y)
- i <- i + 1
- setTxtProgressBar(pb, i)
+ #i <- i + 1
+ #setTxtProgressBar(pb, i)
+ dbDisconnect(db)
+ dev.off()
+ cumulative_data
}
-close(pb)
+#close(pb)
-print(cumulative_data)
-print(summary(cumulative_data))
+rmse_table <- cumulative_data[c("timestamp","x_rmse","y_rmse","z_rmse","speed")]
+write.table(rmse_table, "build/daily-rmse", row.names=FALSE, quote=FALSE)
+#print(cumulative_data)
+#print(summary(cumulative_data))
model <- nls(x_sd ~ a*speed+b,
data=cumulative_data,
control=list(warnOnly=TRUE),
start=list(a=1,b=0))
coefs <- summary(model)$coefficients
-print(coefs)
+#print(coefs)
x_sd_est <- sapply(cumulative_data$speed, function (x) { coefs[[1]]*x + coefs[[2]] })
plot(cumulative_data$speed, cumulative_data$x_sd)
lines(cumulative_data$speed, x_sd_est, col='red')
@@ -470,6 +505,7 @@ lines(cumulative_data$abs_x_mode_acf_x_c, bc_est, col='red')
+
#plot(cumulative_data$x_mode, cumulative_data$c)
#plot(cumulative_data$x_mode, cumulative_data$acf_x_b**cumulative_data$c)
@@ -491,4 +527,3 @@ with (cumulative_data, {
plot(direction, direction(x_b**acf_x_c, y_b**acf_x_c))
})
-dbDisconnect(mydb)
diff --git a/R/calibrate.R b/R/calibrate.R
@@ -1,36 +0,0 @@
-# https://stackoverflow.com/questions/13023274/how-to-do-printf-in-r
-printf <- function(...) cat(sprintf(...))
-
-reference_value <- 2.045 # measured with commercial anemometer
-
-calc_origin <- function (filename="anemometer-no-wind") {
- data <- read.table(filename, sep="", header=FALSE)
- names(data) <- c("t", "x", "y", "z", "sx", "sy", "sz")
- origin = list(x=mean(data$x), y=mean(data$y), z=mean(data$z))
- origin
-}
-
-my_sqrt <- function (x) { sign(x) * sqrt(abs(x)) }
-
-calc_coefficient <- function (filename, axis) {
- origin <- calc_origin(sprintf("%s-origin",filename))
- data <- read.table(filename, sep="", header=FALSE)
- names(data) <- c("t", "x", "y", "z", "sx", "sy", "sz")
- sign <- data[1,sprintf("s%s",axis)]
- sign*my_sqrt(mean(data[[axis]])-origin[[axis]])/reference_value
-}
-
-coefficient_x <- calc_coefficient("build/x", "x")
-coefficient_x_reverse <- calc_coefficient("build/x-reverse", "x")
-coefficient_y <- calc_coefficient("build/y", "y")
-coefficient_y_reverse <- calc_coefficient("build/y-reverse", "y")
-coefficient_z <- calc_coefficient("build/z", "z")
-coefficient_z_reverse <- calc_coefficient("build/z-reverse", "z")
-printf("-1 m/s is x %f y %f z %f\n", coefficient_x, coefficient_y, coefficient_z)
-printf("1 m/s is x %f y %f z %f\n", coefficient_x_reverse, coefficient_y_reverse, coefficient_z_reverse)
-
-weight_0 = (calc_origin("build/z-fan")$z - calc_origin("build/z-fan-origin")$z)
-weight_1 = (calc_origin("build/z-fan-and-wind")$z - calc_origin("build/z-fan-and-wind-origin")$z)
-print(weight_0)
-print(weight_1)
-print(my_sqrt(weight_1-weight_0)/reference_value)
diff --git a/R/common.R b/R/common.R
@@ -0,0 +1,56 @@
+# https://stackoverflow.com/questions/13023274/how-to-do-printf-in-r
+printf <- function(...) cat(sprintf(...))
+
+reference_value <- 2.045 # measured with commercial anemometer
+
+calc_origin <- function (filename="anemometer-no-wind") {
+ data <- read.table(filename, sep="", header=FALSE)
+ names(data) <- c("t", "x", "y", "z", "sx", "sy", "sz")
+ origin = list(x=mean(data$x), y=mean(data$y), z=mean(data$z))
+ origin
+}
+
+abs_sqrt <- function (x) { sign(x) * sqrt(abs(x)) }
+
+calc_coefficient <- function (filename, axis) {
+ origin <- calc_origin(sprintf("%s-origin",filename))
+ data <- read.table(filename, sep="", header=FALSE)
+ names(data) <- c("t", "x", "y", "z", "sx", "sy", "sz")
+ sign <- data[1,sprintf("s%s",axis)]
+ sign*abs_sqrt(mean(data[[axis]])-origin[[axis]])/reference_value
+}
+
+calibration_coefficients <- function(debug=FALSE) {
+ # *1 --- negative, *2 --- positive, *reverse --- positive
+ cx1 <- calc_coefficient("build/x", "x")
+ cx2 <- calc_coefficient("build/x-reverse", "x")
+ cy1 <- calc_coefficient("build/y", "y")
+ cy2 <- calc_coefficient("build/y-reverse", "y")
+ cz1 <- calc_coefficient("build/z", "z")
+ cz2 <- calc_coefficient("build/z-reverse", "z")
+ if (debug) {
+ printf("-1 m/s is x %f y %f z %f\n", cx1, cy1, cz1)
+ printf("1 m/s is x %f y %f z %f\n", cx2, cy2, cz2)
+ }
+ list(x1=cx1, x2=cx2, y1=cy1, y2=cy2, z1=cz1, z2=cz2)
+ #weight_0 = (calc_origin("build/z-fan")$z - calc_origin("build/z-fan-origin")$z)
+ #weight_1 = (calc_origin("build/z-fan-and-wind")$z - calc_origin("build/z-fan-and-wind-origin")$z)
+ #print(weight_0)
+ #print(weight_1)
+ #print(abs_sqrt(weight_1-weight_0)/reference_value)
+}
+
+# convert load cell anemometer sensor value to wind speed in m/s
+sample_to_speed <- function(x, c1, c2) {
+ # remove linear trend
+ t <- c(1:length(x))
+ reg <- lm(x~t)
+ x <- x - reg$fitted.values
+ #x <- x - mean(x) # convert from unsigned to signed sensor values
+ x <- sign(x)*sqrt(abs(x)) # convert from force to velocity
+ # scale sensor values to wind speed using calibration coefficients
+ x[x<0] = x[x<0] / -c1
+ x[x>0] = x[x>0] / c2
+ x
+}
+
diff --git a/R/daily-stats.R b/R/daily-stats.R
@@ -1,19 +1,6 @@
library(DBI)
-source("R/calibrate.R")
-
-sample_to_speed <- function(x, coef, coef_reverse) {
- x <- x - mean(x) # convert from unsigned to signed sensor values
- x <- sign(x)*sqrt(abs(x)) # convert from force to velocity
- # scale sensor values to wind speed using calibration coefficients
- x[x<0] = x[x<0] / -coef
- x[x>0] = x[x>0] / coef_reverse
- # remove linear trend
- t <- c(1:length(x))
- reg <- lm(x~t)
- x <- x - reg$fitted.values
- x
-}
+source("R/common.R")
db <- dbConnect(RSQLite::SQLite(), "samples/load-cell.sqlite3")
RSQLite::initExtension(db)
@@ -30,13 +17,14 @@ i <- 0
result <- data.frame(interval=numeric(0),timestamp=numeric(0),
x_sd=numeric(0),y_sd=numeric(0),z_sd=numeric(0),
speed=numeric(0))
+c <- calibration_coefficients()
unique_intervals <- unique(velocity$interval)
pb <- txtProgressBar(min = 0, max = length(unique_intervals), style = 3)
for (interval in unique_intervals) {
v <- (velocity[velocity$interval==interval, ])
- v$x <- sample_to_speed(v$x, coefficient_x, coefficient_x_reverse)
- v$y <- sample_to_speed(v$y, coefficient_y, coefficient_y_reverse)
- v$z <- sample_to_speed(v$z, coefficient_z, coefficient_z_reverse)
+ v$x <- sample_to_speed(v$x, c$x1, c$x2)
+ v$y <- sample_to_speed(v$y, c$y1, c$y2)
+ v$z <- sample_to_speed(v$z, c$z1, c$z2)
new_row = nrow(result)+1
result[new_row, "interval"] <- interval
result[new_row, "timestamp"] <- min(v$timestamp)
@@ -47,6 +35,7 @@ for (interval in unique_intervals) {
i <- i + 1
setTxtProgressBar(pb, i)
}
+close(pb)
write.table(result, "build/daily-stats", row.names=FALSE, quote=FALSE)
dbDisconnect(db)
diff --git a/R/minmax.R b/R/minmax.R
@@ -0,0 +1,17 @@
+library(DBI)
+
+data <- read.table("build/daily-rmse", sep="", header=TRUE)
+z_rmse_sorted <- sort(data$z_rmse, decreasing=TRUE)
+print("MAX RMSE")
+print(data[data$x_rmse == max(data$x_rmse), ])
+print(data[data$y_rmse == max(data$y_rmse), ])
+print(data[data$z_rmse == z_rmse_sorted[[1]], ])
+print(data[data$z_rmse == z_rmse_sorted[[2]], ])
+print("MIN RMSE")
+print(data[data$x_rmse == min(data$x_rmse), ])
+print(data[data$y_rmse == min(data$y_rmse), ])
+print(data[data$z_rmse == min(data$z_rmse), ])
+print("AVG RMSE")
+print(summary(data))
+print(cor(data))
+#print(data)
diff --git a/gnuplot/daily-rmse.gnuplot b/gnuplot/daily-rmse.gnuplot
@@ -0,0 +1,31 @@
+load 'gnuplot/paired.pal'
+set terminal svg size 500,250 font 'Free Serif, 10' enhanced dynamic round
+set xtics nomirror out rotate 90 offset 0,0.25
+set ytics nomirror out offset 0.5,0
+set border 1+2
+set key top center outside maxrows 1
+
+stats 'build/daily-rmse' using 2
+
+# time axis
+set xdata time
+set timefmt "%s"
+set format x "%Y-%m-%d"
+set xtics 0,60*60*24
+set mxtics 1
+localtime(t) = t+3*60*60
+
+set output 'build/daily-rmse.svg'
+
+do for [timestamp in "1614546000 1614978000 1615496400 1616187600 1616274000 1616360400 1616533200"] {
+ t = int(timestamp)
+ set obj t rectangle behind from first localtime(t), first 0 to first localtime(t)+24*60*60, graph 1 back fillstyle solid 1.0 fillcolor '#f0f0c0' lw 0
+}
+
+plot \
+ 'build/daily-rmse' using (localtime($1)):2 with lines lc '#c04040' title 'RMSE_x',\
+ 'build/daily-rmse' using (localtime($1)):3 with lines lc '#40c040' title 'RMSE_y',\
+ 'build/daily-rmse' using (localtime($1)):4 with lines lc '#4040c0' title 'RMSE_z',\
+ STATS_mean with lines lc '#808080' dashtype 2 notitle
+
+
diff --git a/gnuplot/daily-stats.gnuplot b/gnuplot/daily-stats.gnuplot
@@ -1,5 +1,5 @@
load 'gnuplot/paired.pal'
-set terminal svg size 500,400 font 'Free Serif, 10' enhanced dynamic
+set terminal svg size 500,400 font 'Free Serif, 10' enhanced dynamic round
set xtics nomirror out rotate 90 offset 0,0.25
set ytics nomirror out offset 0.5,0
set border 1+2
diff --git a/gnuplot/velocity-dist-weibull2.gnuplot b/gnuplot/velocity-dist-weibull2.gnuplot
@@ -0,0 +1,107 @@
+load 'gnuplot/paired.pal'
+set terminal svg size 500,500/3*2 font 'Free Serif, 10' enhanced round dynamic
+set xtics nomirror out offset 0,0.25
+set ytics nomirror out offset 0.5,0
+set border 1+2 back
+#set key top center outside Left reverse
+unset key
+
+set output 'build/gnuplot/velocity-dist-weibull2.svg'
+
+# 258 1615914000 0.08700399 0.06517805 0.06370091 1.660843
+# timestamp x_rmse y_rmse z_rmse speed
+# 295 1616187600 0.03527403 0.08648245 0.03217789 1.637611
+# timestamp x_rmse y_rmse z_rmse speed
+# 399 1616958000 0.02882994 0.03539276 0.1744659 0.8539361
+# timestamp x_rmse y_rmse z_rmse speed
+# 292 1616166000 0.06300805 0.02371221 0.094087 2.052675
+# [1] "MIN RMSE"
+# timestamp x_rmse y_rmse z_rmse speed
+# 322 1616382000 0.008555768 0.009697918 0.01404791 3.981209
+# timestamp x_rmse y_rmse z_rmse speed
+# 115 1614877200 0.01305041 0.009301931 0.02022811 4.356801
+# timestamp x_rmse y_rmse z_rmse speed
+# 39 1614330000 0.01867408 0.01648803 0.01133067 4.552804
+
+
+set xlabel 'υ_x' offset 0,1.0
+set yrange [0:1.2]
+set rmargin 0
+set tmargin 1
+set multiplot layout 2,3
+
+positive(x) = x<0 ? (1/0) : x
+negative(x) = x>0 ? (1/0) : x
+
+set lmargin 5
+set title offset 0,-1
+set title sprintf('RMSE_x = %f', 0.08700399)
+set xrange [-4:4]
+plot \
+'build/velocity/x/weibull2/1615914000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_x',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_x'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+#set key top center outside Left reverse
+set lmargin 5
+#set xrange [-3:3]
+set xlabel 'υ_y'
+set title sprintf('RMSE_y = %f', 0.08648245)
+set xrange [-3:5]
+plot \
+'build/velocity/y/weibull2/1616187600' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_y',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_y'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+#unset key
+
+set rmargin 2
+set xrange [-4:4]
+#set xtics -2,1,2
+set xlabel 'υ_z'
+set title sprintf('RMSE_z = %f', 0.1744659)
+plot \
+'build/velocity/z/weibull2/1616958000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_z',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_z',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_z'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+###################################
+# second row
+###################################
+set bmargin 3
+set xrange [-10:10]
+set yrange [0:0.4]
+set ytics 0,0.1
+#set xtics -10,1
+
+set lmargin 5
+set title sprintf('RMSE_x = %f', 0.008555768)
+plot \
+'build/velocity/x/weibull2/1616382000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_x',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_x'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+set lmargin 5
+#set xrange [-3:3]
+set xlabel 'υ_y'
+set title sprintf('RMSE_y = %f', 0.009301931)
+#set xrange [-3:5]
+plot \
+'build/velocity/y/weibull2/1614877200' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_y',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_y'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+set rmargin 2
+#set xtics -2,1,2
+set xlabel 'υ_z'
+set title sprintf('RMSE_z = %f', 0.01133067)
+plot \
+'build/velocity/z/weibull2/1614330000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_z',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_z',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_z'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
diff --git a/gnuplot/velocity-dist.gnuplot b/gnuplot/velocity-dist.gnuplot
@@ -0,0 +1,107 @@
+load 'gnuplot/paired.pal'
+set terminal svg size 500,500/3*2 font 'Free Serif, 10' enhanced round dynamic
+set xtics nomirror out offset 0,0.25
+set ytics nomirror out offset 0.5,0
+set border 1+2 back
+#set key top center outside Left reverse
+unset key
+
+set output 'build/gnuplot/velocity-dist.svg'
+
+# 258 1615914000 0.08700399 0.06517805 0.06370091 1.660843
+# timestamp x_rmse y_rmse z_rmse speed
+# 295 1616187600 0.03527403 0.08648245 0.03217789 1.637611
+# timestamp x_rmse y_rmse z_rmse speed
+# 399 1616958000 0.02882994 0.03539276 0.1744659 0.8539361
+# timestamp x_rmse y_rmse z_rmse speed
+# 292 1616166000 0.06300805 0.02371221 0.094087 2.052675
+# [1] "MIN RMSE"
+# timestamp x_rmse y_rmse z_rmse speed
+# 322 1616382000 0.008555768 0.009697918 0.01404791 3.981209
+# timestamp x_rmse y_rmse z_rmse speed
+# 115 1614877200 0.01305041 0.009301931 0.02022811 4.356801
+# timestamp x_rmse y_rmse z_rmse speed
+# 39 1614330000 0.01867408 0.01648803 0.01133067 4.552804
+
+
+set xlabel 'υ_x' offset 0,1.0
+set yrange [0:1.2]
+set rmargin 0
+set tmargin 1
+set multiplot layout 2,3
+
+positive(x) = x<0 ? (1/0) : x
+negative(x) = x>0 ? (1/0) : x
+
+set lmargin 5
+set title offset 0,-1
+set title sprintf('RMSE_x = %f', 0.08700399)
+set xrange [-4:4]
+plot \
+'build/velocity/x/weibull/1615914000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_x',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_x'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+#set key top center outside Left reverse
+set lmargin 5
+#set xrange [-3:3]
+set xlabel 'υ_y'
+set title sprintf('RMSE_y = %f', 0.08648245)
+set xrange [-3:5]
+plot \
+'build/velocity/y/weibull/1616187600' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_y',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_y'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+#unset key
+
+set rmargin 2
+set xrange [-4:4]
+#set xtics -2,1,2
+set xlabel 'υ_z'
+set title sprintf('RMSE_z = %f', 0.1744659)
+plot \
+'build/velocity/z/weibull/1616958000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_z',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_z',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_z'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+###################################
+# second row
+###################################
+set bmargin 3
+set xrange [-10:10]
+set yrange [0:0.4]
+set ytics 0,0.1
+#set xtics -10,1
+
+set lmargin 5
+set title sprintf('RMSE_x = %f', 0.008555768)
+plot \
+'build/velocity/x/weibull/1616382000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_x',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_x'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+set lmargin 5
+#set xrange [-3:3]
+set xlabel 'υ_y'
+set title sprintf('RMSE_y = %f', 0.009301931)
+#set xrange [-3:5]
+plot \
+'build/velocity/y/weibull/1614877200' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_y',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_y'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
+set rmargin 2
+#set xtics -2,1,2
+set xlabel 'υ_z'
+set title sprintf('RMSE_z = %f', 0.01133067)
+plot \
+'build/velocity/z/weibull/1614330000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_z',\
+'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_z',\
+'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_z'
+#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+
diff --git a/guix/manifest.scm b/guix/manifest.scm
@@ -4,6 +4,8 @@
(@ (gnu packages statistics) r-rsqlite)
(@ (gnu packages statistics) r-sn)
(@ (gnu packages statistics) r-plotrix)
+ (@ (gnu packages statistics) r-foreach)
+ (@ (gnu packages statistics) r-doparallel)
(@ (gnu packages cran) r-fitdistrplus)
(@ (gnu packages statistics) r)
(@ (gnu packages maths) gnuplot)
diff --git a/main.tex b/main.tex
@@ -237,15 +237,15 @@ listing~\ref{lst-sample-to-speed}.
\begin{lstlisting}[label={lst-sample-to-speed},caption={The code that transforms raw load cell sensor values into wind speed projection to the corresponding axis.}]
sampleToSpeed <- function(x, c1, c2) {
- x <- x - mean(x) # convert from unsigned to signed sensor values
- x <- sign(x)*sqrt(abs(x)) # convert from force to velocity
- # scale sensor values to wind speed using calibration coefficients
- x[x<0] = x[x<0] / c1
- x[x>0] = x[x>0] / c2
# remove linear trend
t <- c(1:length(x))
reg <- lm(x~t)
x - reg$fitted.values
+ x <- sign(x)*sqrt(abs(x)) # convert from force to velocity
+ # scale sensor values to wind speed using calibration coefficients
+ x[x<0] = x[x<0] / c1
+ x[x>0] = x[x>0] / c2
+ x
}
\end{lstlisting}
@@ -315,41 +315,37 @@ statistics that is not distorted by the turbulence.
\subsection{Verification of wind velocity distribution}
-The data collected with three-axis anemometer showed that for each two-hour
-interval the shape of the wind velocity distribution for each axis can be
-approximated by~\eqref{eq-velocity-distribution}. If the wind blows in positive
-\(x\) direction, then the right part of the distribution is steep and if the
-wind blows in negative \(x\) direction, then the left part of the distribution
-is steep.
-
-We failed to gather the data where the distribution is unimodal: in most
-intervals the distribution is bimodal. However, the second peak is often small
-enough to have small effect on least-squares fitting. Distributions for each
-axis are presented in figure~\ref{fig-velocity-distributions} and their
-parameters in table~\ref{tab-velocity-distributions}.
+The wind speed data collected with three-axis anemometer was approimated by
+Weibull distribution using least-squares fitting. Negative and positive wind
+speed projections to each axis both have this distribution, but with different
+parameters. Most of the data intervals contain only one mean wind direction, which
+means that one of the distributions is for incident wind flow on the arm of the
+anemometer and another one is for the turbulent flow that forms behind the arm.
+For \(z\) axis both left and right distributions have similar shapes, for \(x\)
+and \(y\) axes the distribution for incident flow is taller than the
+distribution for turbulent flow. Distributions for each axis are presented in
+figure~\ref{fig-velocity-distributions}.
+
+RMSE of the approximation has negative correlation with wind speed: the larger
+the wind speed, the smaller the error and vice versa. Larger error for low wind
+speeds is caused by larger skewness and kurtosis (see the first row of
+figure~\ref{fig-velocity-distributions}). RMSE is reduced
+when~\eqref{eq-velocity-distribution} is used as the distribution for both low
+and high wind speeds. TODO
\begin{figure}
\centering
- TODO
- \caption{Per-axis wind velocity distributions fitted into the model
- described by~\eqref{eq-velocity-distribution}.
+ \includegraphics{build/gnuplot/velocity-dist.eps}
+ \caption{Per-axis wind velocity distributions fitted into Weibull
+ distribution.
Data intervals with the smallest error (the first row),
- data intervals with the largest error (the second row).
+ data intervals with the largest error (the second row). Red line shows estimated
+ probability density of positive wind speed projections, blue line shows estimated
+ probability density of negative wind speed projections and circles denote observed
+ probability density of wind speed projections.
\label{fig-velocity-distributions}}
\end{figure}
-\begin{table}
- \centering
- \begin{tabular}{ll}
- \toprule
- TODO \\
- \bottomrule
- \end{tabular}
- \caption{Model fitting errors for each graph in
- figure~\ref{fig-velocity-distributions}.\label{tab-velocity-distributions}}
-\end{table}
-
-
\subsection{Verification of wind direction distribution}
\cite{carta2008}
