commit aa45203011b608cfc43c042ee366023631245988
parent 8011b596532b33ffd4a95162e2c61c8e20c8f7a3
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Mon, 3 May 2021 16:14:32 +0300
Calculate turbulence.
Diffstat:
4 files changed, 143 insertions(+), 53 deletions(-)
diff --git a/Makefile b/Makefile
@@ -24,6 +24,7 @@ build/main.pdf: build/gnuplot/direction-dist.eps
build/main.pdf: build/gnuplot/acf.eps
build/main.pdf: build/gnuplot/hold-peak.eps
build/main.pdf: build/gnuplot/vtestbed.eps
+build/main.pdf: build/gnuplot/turbulence.eps
#build/main.pdf: build/gnuplot/anemometer.eps
#build/main.pdf: build/gnuplot/anemometer.svg
build/main.pdf:
@@ -68,6 +69,7 @@ build/gnuplot/velocity-xy-dist-max.svg: gnuplot/velocity-xy-dist.gnuplot build/g
build/gnuplot/direction-dist.svg: build/gnuplot/direction_rmse.gnuplot
build/gnuplot/vtestbed.svg: gnuplot/vtestbed.gnuplot
+build/gnuplot/turbulence.svg: gnuplot/turbulence.gnuplot
#build/main.zip: build/gnuplot/*.eps main.tex
build/main.zip: main.tex
diff --git a/R/anal.R b/R/anal.R
@@ -9,6 +9,11 @@ registerDoParallel(cores=detectCores())
source("R/common.R")
source("R/distributions.R")
+# https://stackoverflow.com/questions/4954507/calculate-the-area-under-a-curve
+area_under_curve <- function (x,y) {
+ sum(diff(x) * (head(y,-1)+tail(y,-1)))/2
+}
+
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;
@@ -36,8 +41,12 @@ fit_velocity_pdf_bilateral <- function (velocity, axis="x", timestamp=0, dist="w
write("\n",file=filename,append=TRUE)
write.table(data.frame(x=d$x,density=d$y), filename,
row.names=FALSE, quote=FALSE, append=TRUE)
+ area_left <- area_under_curve(dist_left$x, dist_left$actual_density)
+ area_right <- area_under_curve(dist_right$x, dist_right$actual_density)
list(coefficients=dist_left$coefficients,
- rmse=rmse(estimated_density, actual_density))
+ rmse=rmse(estimated_density, actual_density),
+ area_min=min(area_left,area_right),
+ area_max=max(area_left,area_right))
}
fit_velocity_acf_bilateral <- function (velocity, axis="x", timestamp=0, plot=TRUE) {
@@ -117,10 +126,14 @@ fit_direction_bilateral <- function (direction, plot=TRUE, timestamp=0, dist="vo
time_delta = 60*60*2
#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_start = 1616936400
-time_end = time_start
+
+# all of the data
+time_start = time_string_to_timestamp("2021-02-21T00:00")
+time_end = time_string_to_timestamp("2021-03-31T23:00")
+
+# timestamp for ACF validation
+#time_start = 1616936400
+#time_end = time_start
print(time_start)
print(time_end)
@@ -150,7 +163,10 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
acf_y_a=numeric(0),acf_y_b=numeric(0),acf_y_c=numeric(0),
acf_z_a=numeric(0),acf_z_b=numeric(0),acf_z_c=numeric(0),
acf_v_a=numeric(0),acf_v_b=numeric(0),acf_v_c=numeric(0),
- acf_d_a=numeric(0),acf_v_b=numeric(0),acf_v_c=numeric(0))
+ acf_d_a=numeric(0),acf_v_b=numeric(0),acf_v_c=numeric(0),
+ x_area_min=numeric(0),x_area_max=numeric(0),
+ y_area_min=numeric(0),y_area_max=numeric(0),
+ z_area_min=numeric(0),z_area_max=numeric(0))
velocity <- select_samples(timestamp, time_delta)
if (nrow(velocity) == 0) { return(NULL) }
filename <- sprintf("build/rplot-%d.pdf", timestamp)
@@ -166,6 +182,40 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
cumulative_data[new_row, "z_mode"] = z_mode
cumulative_data[new_row, "speed"] = mean(velocity$speed)
cumulative_data[new_row, "direction"] = direction(x_mode, y_mode)
+
+ max_abs_mean <- function (x) {
+ #max(abs(mean(x[x<=0])), abs(mean(x[x>=0])))
+ mean_1 <- abs(mean(x[x<=0]))
+ mean_2 <- abs(mean(x[x>=0]))
+ if (mean_1 > mean_2) {
+ return(mean(x[x<=0]))
+ } else {
+ return(mean(x[x>=0]))
+ }
+ }
+
+ max_abs_var <- function (x) {
+ mean_1 <- abs(mean(x[x<=0]))
+ mean_2 <- abs(mean(x[x>=0]))
+ if (mean_1 > mean_2) {
+ return(var(x[x<=0]))
+ } else {
+ return(var(x[x>=0]))
+ }
+ }
+
+ cumulative_data[new_row, "abs_x_mode"] = abs(x_mode)
+ cumulative_data[new_row, "x_sd"] = sd(velocity$x)
+ cumulative_data[new_row, "y_sd"] = sd(velocity$y)
+ cumulative_data[new_row, "z_sd"] = sd(velocity$z)
+ cumulative_data[new_row, "x_mean"] = max_abs_mean(velocity$x)
+ cumulative_data[new_row, "y_mean"] = max_abs_mean(velocity$y)
+ cumulative_data[new_row, "z_mean"] = max_abs_mean(velocity$z)
+ cumulative_data[new_row, "x_var"] = max_abs_var(velocity$x)
+ cumulative_data[new_row, "y_var"] = max_abs_var(velocity$y)
+ cumulative_data[new_row, "z_var"] = max_abs_var(velocity$z)
+ cumulative_data[new_row, "v_mean"] = mean(velocity$speed)
+ cumulative_data[new_row, "d_mean"] = mean(velocity$direction)
#cumulative_data[new_row, "direction"] = weighted.mean(velocity$direction, speed(velocity$x,velocity$y))
#fit <- fitdist(velocity$x, "sn", start=list(omega=1,alpha=-1))
#fit <- fitdist(velocity$x, "weibull")
@@ -183,6 +233,8 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
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_area_min"] = x_model$area_min
+ cumulative_data[new_row, "x_area_max"] = x_model$area_max
#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]]
@@ -192,6 +244,8 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
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_area_min"] = y_model$area_min
+ cumulative_data[new_row, "y_area_max"] = y_model$area_max
#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]]
@@ -199,6 +253,8 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
# Velocity Z
z_model <- fit_velocity_pdf_bilateral(velocity$z, axis="z", timestamp=timestamp, dist="weibull")
cumulative_data[new_row, "z_rmse"] = z_model$rmse
+ cumulative_data[new_row, "z_area_min"] = z_model$area_min
+ cumulative_data[new_row, "z_area_max"] = z_model$area_max
# ACF X
x_model <- fit_velocity_acf_bilateral(velocity$x, plot=TRUE, axis="x", timestamp=timestamp)
cumulative_data[new_row, "acf_x_a"] = x_model$coefficients[[1]]
@@ -255,40 +311,6 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
polar.plot(direction_hist$direction,rad2deg(direction_hist$angle),
start=45,clockwise=TRUE,rp.type="polygon")
par(mar=old_mar)
-
- max_abs_mean <- function (x) {
- #max(abs(mean(x[x<=0])), abs(mean(x[x>=0])))
- mean_1 <- abs(mean(x[x<=0]))
- mean_2 <- abs(mean(x[x>=0]))
- if (mean_1 > mean_2) {
- return(mean(x[x<=0]))
- } else {
- return(mean(x[x>=0]))
- }
- }
-
- max_abs_var <- function (x) {
- mean_1 <- abs(mean(x[x<=0]))
- mean_2 <- abs(mean(x[x>=0]))
- if (mean_1 > mean_2) {
- return(var(x[x<=0]))
- } else {
- return(var(x[x>=0]))
- }
- }
-
- cumulative_data[new_row, "abs_x_mode"] = abs(x_mode)
- cumulative_data[new_row, "x_sd"] = sd(velocity$x)
- cumulative_data[new_row, "y_sd"] = sd(velocity$y)
- cumulative_data[new_row, "z_sd"] = sd(velocity$z)
- cumulative_data[new_row, "x_mean"] = max_abs_mean(velocity$x)
- cumulative_data[new_row, "y_mean"] = max_abs_mean(velocity$y)
- cumulative_data[new_row, "z_mean"] = max_abs_mean(velocity$z)
- cumulative_data[new_row, "x_var"] = max_abs_var(velocity$x)
- cumulative_data[new_row, "y_var"] = max_abs_var(velocity$y)
- cumulative_data[new_row, "z_var"] = max_abs_var(velocity$z)
- cumulative_data[new_row, "v_mean"] = mean(velocity$speed)
- cumulative_data[new_row, "d_mean"] = mean(velocity$direction)
#i <- i + 1
#setTxtProgressBar(pb, i)
dev.off()
@@ -317,6 +339,31 @@ lines(cumulative_data$speed, x_sd_est, col='red')
#plot(log(cumulative_data$acf_x_b), cumulative_data$x_sd)
#plot(cumulative_data$acf_x_c, cumulative_data$x_sd)
+# turbulence coefficient
+#plot(hist(cumulative_data$x_area_min))
+#plot(hist(cumulative_data$x_area_max))
+turbulence_x = cumulative_data$x_area_min / cumulative_data$x_area_max
+turbulence_y = cumulative_data$y_area_min / cumulative_data$y_area_max
+turbulence_z = cumulative_data$z_area_min / cumulative_data$z_area_max
+plot(hist(turbulence_x, breaks=50))
+plot(abs(cumulative_data$x_mean),turbulence_x)
+reg <- lm(turbulence_x~cumulative_data$x_mean)
+lines(cumulative_data$x_mean, reg$fitted.values, col='red')
+print(reg)
+print(mean(turbulence_x))
+#plot(loess(turbulence_x ~ abs(x_mean), cumulative_data))
+scatter.smooth(abs(cumulative_data$x_mean), turbulence_x)
+path <- file.path('build', 'turbulence')
+write.table(data.frame(velocity_x=abs(cumulative_data$x_mean),
+ velocity_y=abs(cumulative_data$y_mean),
+ velocity_z=abs(cumulative_data$z_mean),
+ turbulence_x=turbulence_x,
+ turbulence_y=turbulence_y,
+ turbulence_z=turbulence_z),
+ path, row.names=FALSE, quote=FALSE)
+#lo <- loess.smooth(cumulative_data$x_mean, cumulative_data$turbulence_x, span=0.5)
+#lines(lo$x, lo$y)
+
# remove outliers
cumulative_data <- cumulative_data[cumulative_data$acf_x_b<2 & cumulative_data$acf_y_b<2 & cumulative_data$acf_z_b<2,]
diff --git a/R/vtestbed.R b/R/vtestbed.R
@@ -10,9 +10,9 @@ time_start = 1616936400
time_end = time_start
time_delta = 60*60*2
-#data <- read.table("samples/vtestbed/statistics.log", sep=",", header=TRUE)
+data <- read.table("samples/vtestbed/statistics.log", sep=",", header=TRUE)
#data <- read.table("samples/vtestbed/statistics-2021-04-28.log", sep=",", header=TRUE)
-data <- read.table("samples/vtestbed/statistics-2021-04-28_21-42-37.log", sep=",", header=TRUE)
+#data <- read.table("samples/vtestbed/statistics-2021-04-28_21-42-37.log", sep=",", header=TRUE)
#data <- read.table("samples/vtestbed/wind_data", sep=",", header=TRUE)
#timestamps <- seq(1,floor(max(data$time)))
##print(timestamps)
@@ -113,20 +113,21 @@ banner("Anemometer")
data_real <- select_samples(time_start, time_delta)
dist_speed_real <- do_fit_velocity(data_real$speed, suffix="anem")
dist_direction_real <- do_fit_direction(data_real$direction[data_real$direction<=0], suffix="anem")
-model_acf_x_real <- do_fit_acf(data_real$x[data_real$x<=0], axis="X", suffix="anem", delevel=TRUE)
-model_acf_y_real <- do_fit_acf(data_real$y[data_real$y<=0], axis="Y", suffix="anem", delevel=TRUE)
-model_acf_z_real <- do_fit_acf(data_real$z[data_real$z<=0], axis="Z", suffix="anem", delevel=TRUE)
+model_acf_x_real <- do_fit_acf(data_real$x[data_real$x<=0], axis="X", suffix="anem", delevel=FALSE)
+model_acf_y_real <- do_fit_acf(data_real$y[data_real$y<=0], axis="Y", suffix="anem", delevel=FALSE)
+model_acf_z_real <- do_fit_acf(data_real$z[data_real$z<=0], axis="Z", suffix="anem", delevel=FALSE)
# https://stackoverflow.com/questions/27455948/in-r-whats-the-simplest-way-to-scale-a-vector-to-a-unit-vector
normalise <- function(x) {
- x <- as.vector(x)
- s <- max(abs(x))
- if (s == 0) {
- ret <- x
- } else {
- ret <- x/s
- }
- ret
+ #x <- as.vector(x)
+ #s <- max(abs(x))
+ #if (s == 0) {
+ # ret <- x
+ #} else {
+ # ret <- x/s
+ #}
+ #ret
+ x
}
banner("Comparison")
diff --git a/gnuplot/turbulence.gnuplot b/gnuplot/turbulence.gnuplot
@@ -0,0 +1,40 @@
+set terminal svg size 450,450/3 font 'Free Serif, 10' enhanced round dynamic
+set xtics nomirror out offset 0,0.5
+set ytics nomirror out offset 0.5,0
+set border 1+2 back
+set xlabel offset 0,1.25
+set ylabel offset 2.5,0
+set title offset 0,-1
+set tmargin 1
+set lmargin 6
+set rmargin 1
+unset key
+
+set output 'build/gnuplot/turbulence.svg'
+
+set multiplot layout 1,3
+
+set xrange [0:6]
+set yrange [0.3:*]
+set xtics 0,1
+set ytics 0.4,0.2
+
+n = 3
+
+set xlabel 'Velocity X, m/s'
+set ylabel 'Turbulence coefficient X'
+plot \
+'< sort -nk1 build/turbulence' using (abs($1)):4 every n with points pt 6 lc '#404040' notitle,\
+'' using (abs($1)):4 smooth bezier with lines lc '#c04040' notitle
+
+set xlabel 'Velocity Y, m/s'
+set ylabel 'Turbulence coefficient Y'
+plot \
+'< sort -nk2 build/turbulence' using (abs($2)):5 every n with points pt 6 lc '#404040' notitle,\
+'' using (abs($2)):5 smooth bezier with lines lc '#c04040' notitle
+
+set xlabel 'Velocity Z, m/s'
+set ylabel 'Turbulence coefficient Z'
+plot \
+'< sort -nk3 build/turbulence' using (abs($3)):6 every n with points pt 6 lc '#404040' notitle,\
+'' using (abs($3)):6 smooth bezier with lines lc '#c04040' notitle
