vtestbed.R (6741B)
1 source("R/common.R") 2 source("R/distributions.R") 3 4 delevel <- function (x) { 5 x - min(x) 6 } 7 8 # the timestamp against which the simulation output is validated 9 time_start = 1616936400 10 time_end = time_start 11 time_delta = 60*60*2 12 13 data <- read.table("samples/vtestbed/statistics.log", sep=",", header=TRUE) 14 #data <- read.table("samples/vtestbed/statistics-2021-04-28.log", sep=",", header=TRUE) 15 #data <- read.table("samples/vtestbed/statistics-2021-04-28_21-42-37.log", sep=",", header=TRUE) 16 #data <- read.table("samples/vtestbed/wind_data", sep=",", header=TRUE) 17 #timestamps <- seq(1,floor(max(data$time))) 18 ##print(timestamps) 19 ## resample the data to 1 second period 20 #new_data <- data.frame( 21 # time=timestamps, 22 # wind_velocity_x = approx(data$time, data$wind_velocity_x, timestamps)$y, 23 # wind_velocity_y = approx(data$time, data$wind_velocity_y, timestamps)$y, 24 # wind_velocity_z = approx(data$time, data$wind_velocity_z, timestamps)$y 25 #) 26 #data = new_data 27 28 # select every 10th point 29 #data <- data[seq(1,nrow(data),4000),] 30 31 data$time = seq(1,nrow(data)) 32 x_mean = -2.439 33 y_mean = -2.158 34 z_mean = -1.367 35 data$wind_velocity_x = data$wind_velocity_x + x_mean 36 data$wind_velocity_y = data$wind_velocity_y + y_mean 37 data$wind_velocity_z = data$wind_velocity_z + z_mean 38 #data <- data[c("time","wind_velocity_x","wind_velocity_y","wind_velocity_z")] 39 #plot(data$time, data$wind_velocity_x) 40 #plot(data$time, data$wind_velocity_y) 41 #plot(data$time, data$wind_velocity_z) 42 data$speed <- speed(data$wind_velocity_x,data$wind_velocity_y,data$wind_velocity_z) 43 data$direction <- direction(data$wind_velocity_x,data$wind_velocity_y) 44 #hist(data$speed) 45 #hist(data$direction) 46 47 48 do_fit_velocity <- function (velocity, suffix="no-suffix") { 49 velocity_hist <- normalize_hist(hist(velocity, plot=FALSE, breaks=100)) 50 x = velocity_hist$x; 51 v = velocity_hist$density; 52 d = density(velocity,adjust=0.5) 53 dist <- fit_velocity_pdf(x, v, density=d$y, sign=1, axis="v", dist="weibull") 54 plot(x, dist$actual_density, xlab='Wind speed', ylab='PDF') 55 lines(d, col='green', lwd=2) 56 lines(dist$x, dist$estimated_density, col='red', lwd=2) 57 path <- file.path('build', 'vtestbed', suffix) 58 make_directory(path) 59 filename <- file.path(path, "velocity") 60 write.table(data.frame(x=dist$x,actual=dist$actual_density,estimated=dist$estimated_density), 61 filename, row.names=FALSE, quote=FALSE) 62 dist 63 } 64 65 do_fit_direction <- function (direction, suffix="no-suffix") { 66 hist <- normalize_hist(hist(direction, plot=FALSE, breaks=1000)) 67 x = hist$x; 68 v = hist$density; 69 d = density(direction,adjust=0.5) 70 dist <- fit_direction(x, v, plot=FALSE) 71 plot(x, dist$actual_density, xlab='Direction', ylab='PDF') 72 lines(d, col='green', lwd=2) 73 lines(dist$x, dist$estimated_density, col='red', lwd=2) 74 path <- file.path('build', 'vtestbed', suffix) 75 make_directory(path) 76 filename <- file.path(path, "direction") 77 write.table(data.frame(x=dist$x,actual=dist$actual_density,estimated=dist$estimated_density), 78 filename, row.names=FALSE, quote=FALSE) 79 dist 80 } 81 82 do_fit_acf <- function (velocity, axis="x", suffix="no-suffix", delevel=FALSE) { 83 v = velocity 84 model <- fit_velocity_acf(v, axis=axis, plot=FALSE, write=FALSE); 85 if (delevel) { 86 model$actual = delevel(model$actual) 87 model$estimated = delevel(model$estimated) 88 } 89 plot(model$lag, model$actual, xlab='Lag', ylab=paste('ACF ', axis, sep="")) 90 lines(model$lag, model$estimated, col='red', lwd=2) 91 path <- file.path('build', 'vtestbed', suffix) 92 make_directory(path) 93 filename <- file.path(path, sprintf("acf-%s", axis)) 94 write.table(data.frame(lag=model$lag,actual=model$actual,estimated=model$estimated), 95 filename, row.names=FALSE, quote=FALSE) 96 model 97 } 98 99 # https://stackoverflow.com/questions/19918985/r-plot-only-text 100 banner <- function (str) { 101 plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n') 102 text(x = 0.5, y = 0.5, str, cex = 1.6, col = "black") 103 } 104 105 banner("Virtual testbed") 106 dist_speed <- do_fit_velocity(data$speed, suffix="vtb") 107 dist_direction <- do_fit_direction(data$direction[data$direction<=0], suffix="vtb") 108 model_acf_x <- do_fit_acf(data$wind_velocity_x, axis="X", suffix="vtb") 109 model_acf_y <- do_fit_acf(data$wind_velocity_y, axis="Y", suffix="vtb") 110 model_acf_z <- do_fit_acf(data$wind_velocity_z, axis="Z", suffix="vtb") 111 112 banner("Anemometer") 113 data_real <- select_samples(time_start, time_delta) 114 dist_speed_real <- do_fit_velocity(data_real$speed, suffix="anem") 115 dist_direction_real <- do_fit_direction(data_real$direction[data_real$direction<=0], suffix="anem") 116 model_acf_x_real <- do_fit_acf(data_real$x[data_real$x<=0], axis="X", suffix="anem", delevel=FALSE) 117 model_acf_y_real <- do_fit_acf(data_real$y[data_real$y<=0], axis="Y", suffix="anem", delevel=FALSE) 118 model_acf_z_real <- do_fit_acf(data_real$z[data_real$z<=0], axis="Z", suffix="anem", delevel=FALSE) 119 120 # https://stackoverflow.com/questions/27455948/in-r-whats-the-simplest-way-to-scale-a-vector-to-a-unit-vector 121 normalise <- function(x) { 122 #x <- as.vector(x) 123 #s <- max(abs(x)) 124 #if (s == 0) { 125 # ret <- x 126 #} else { 127 # ret <- x/s 128 #} 129 #ret 130 x 131 } 132 133 banner("Comparison") 134 plot(dist_speed$x, normalise(dist_speed$actual_density), xlab='Wind speed', ylab='PDF') 135 lines(dist_speed$x, normalise(dist_speed$estimated_density), col='red', lwd=2) 136 points(dist_speed_real$x, normalise(dist_speed_real$actual_density), col='blue') 137 lines(dist_speed_real$x, normalise(dist_speed_real$estimated_density), col='blue', lwd=2) 138 139 plot(dist_direction$x, normalise(dist_direction$actual_density), xlab='Wind direction', ylab='PDF') 140 lines(dist_direction$x, normalise(dist_direction$estimated_density), col='red', lwd=2) 141 points(dist_direction_real$x, normalise(dist_direction_real$actual_density), col='blue') 142 lines(dist_direction_real$x, normalise(dist_direction_real$estimated_density), col='blue', lwd=2) 143 144 plot(model_acf_x$lag, model_acf_x$actual, col='red', xlab='Lag', ylab="ACF X") 145 lines(model_acf_x$lag, model_acf_x$estimated, col='red', lwd=2) 146 points(model_acf_x_real$lag, (model_acf_x_real$actual), col='blue') 147 lines(model_acf_x_real$lag, (model_acf_x_real$estimated), col='blue', lwd=2) 148 149 plot(model_acf_y$lag, model_acf_y$actual, col='red', xlab='Lag', ylab="ACF Y") 150 lines(model_acf_y$lag, model_acf_y$estimated, col='red', lwd=2) 151 points(model_acf_y_real$lag, (model_acf_y_real$actual), col='blue') 152 lines(model_acf_y_real$lag, (model_acf_y_real$estimated), col='blue', lwd=2) 153 154 plot(model_acf_z$lag, model_acf_z$actual, col='red', xlab='Lag', ylab="ACF Z") 155 lines(model_acf_z$lag, model_acf_z$estimated, col='red', lwd=2) 156 points(model_acf_z_real$lag, (model_acf_z_real$actual), col='blue') 157 lines(model_acf_z_real$lag, (model_acf_z_real$estimated), col='blue', lwd=2)