iccsa-21-wind

git clone https://git.igankevich.com/iccsa-21-wind.git
Log | Files | Refs

commit 0fda1b187f09aed54ecbbb5c3268e88d462da040
parent 39c7a388abf35e1f90074d2151d6c1d695806f57
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date:   Mon, 26 Apr 2021 16:06:53 +0300

Add comparison graph.

Diffstat:
Makefile | 2++
R/distributions.R | 20+++++++++++++-------
R/vtestbed.R | 79+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------
gnuplot/vtestbed.gnuplot | 66++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
main.bib | 13+++++++++++++
main.tex | 48++++++++++++++++++++++++++++++++++++++----------
6 files changed, 197 insertions(+), 31 deletions(-)

diff --git a/Makefile b/Makefile @@ -23,6 +23,7 @@ build/main.pdf: build/gnuplot/velocity-xy-dist-min.eps 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/anemometer.eps #build/main.pdf: build/gnuplot/anemometer.svg build/main.pdf: @@ -66,6 +67,7 @@ build/gnuplot/velocity-xy-dist-max.svg: gnuplot/velocity-xy-dist.gnuplot build/g gnuplot -d -c gnuplot/velocity-xy-dist.gnuplot max build/gnuplot/direction-dist.svg: build/gnuplot/direction_rmse.gnuplot +build/gnuplot/vtestbed.svg: gnuplot/vtestbed.gnuplot #build/main.zip: build/gnuplot/*.eps main.tex build/main.zip: main.tex diff --git a/R/distributions.R b/R/distributions.R @@ -216,13 +216,19 @@ fit_velocity_acf <- function (velocity, plot=TRUE, axis="x", timestamp=0, write= if (debug) { plot(acf_x$lag, acf_x$acf, xlab='Lag', ylab=paste('DEBUG ACF', axis)) } - model <- nls(acf ~ a*exp(-abs(b*lag)**c), - data=acf_x, - start=list(a=1,b=1,c=1), - control=list(warnOnly=TRUE), - algorithm="port", - lower=c(0,0,1e-6), upper=c(1000,1000,10)) - c <- summary(model)$coefficients + v_mean <- mean(velocity) + c <- tryCatch({ + model <- nls(acf ~ a*exp(-abs(b*lag)**c), + data=acf_x, + start=list(a=1,b=1,c=1), + control=list(warnOnly=TRUE), + algorithm="port", + lower=c(0,0,1e-6), upper=c(1000,1000,10)) + c <- summary(model)$coefficients + },error=function(cond) { + message(cond) + c(1,1,1) + }) acf_est <- c[[1]]*exp(-abs(c[[2]]*acf_x$lag)**c[[3]]) if (plot) { plot(acf_x$lag, acf_x$acf, xlab='Lag', ylab=paste('ACF', axis)) diff --git a/R/vtestbed.R b/R/vtestbed.R @@ -23,7 +23,7 @@ data$direction <- direction(data$wind_velocity_x,data$wind_velocity_y) #hist(data$direction) -do_fit_velocity <- function (velocity) { +do_fit_velocity <- function (velocity, suffix="no-suffix") { velocity_hist <- normalize_hist(hist(velocity, plot=FALSE, breaks=100)) x = velocity_hist$x; v = velocity_hist$density; @@ -32,10 +32,15 @@ do_fit_velocity <- function (velocity) { plot(x, dist$actual_density, xlab='Wind speed', ylab='PDF') lines(d, col='green', lwd=2) lines(dist$x, dist$estimated_density, col='red', lwd=2) + path <- file.path('build', 'vtestbed', suffix) + make_directory(path) + filename <- file.path(path, "velocity") + write.table(data.frame(x=dist$x,actual=dist$actual_density,estimated=dist$estimated_density), + filename, row.names=FALSE, quote=FALSE) dist } -do_fit_direction <- function (direction) { +do_fit_direction <- function (direction, suffix="no-suffix") { hist <- normalize_hist(hist(direction, plot=FALSE, breaks=1000)) x = hist$x; v = hist$density; @@ -44,14 +49,24 @@ do_fit_direction <- function (direction) { plot(x, dist$actual_density, xlab='Direction', ylab='PDF') lines(d, col='green', lwd=2) lines(dist$x, dist$estimated_density, col='red', lwd=2) + path <- file.path('build', 'vtestbed', suffix) + make_directory(path) + filename <- file.path(path, "direction") + write.table(data.frame(x=dist$x,actual=dist$actual_density,estimated=dist$estimated_density), + filename, row.names=FALSE, quote=FALSE) dist } -do_fit_acf <- function (velocity, axis="x") { +do_fit_acf <- function (velocity, axis="x", suffix="no-suffix") { v = velocity model <- fit_velocity_acf(v, axis=axis, plot=FALSE, write=FALSE); plot(model$lag, model$actual, xlab='Lag', ylab=paste('ACF ', axis, sep="")) lines(model$lag, model$estimated, col='red', lwd=2) + path <- file.path('build', 'vtestbed', suffix) + make_directory(path) + filename <- file.path(path, sprintf("acf-%s", axis)) + write.table(data.frame(lag=model$lag,actual=model$actual,estimated=model$estimated), + filename, row.names=FALSE, quote=FALSE) model } @@ -62,18 +77,54 @@ banner <- function (str) { } banner("Virtual testbed") -dist_speed <- do_fit_velocity(data$speed) -dist_direction <- do_fit_direction(data$direction) -#model_acf_x <- do_fit_acf(data$wind_velocity_x, axis="X") -model_acf_y <- do_fit_acf(data$wind_velocity_y, axis="Y") -#model_acf_z <- do_fit_acf(data$wind_velocity_z, axis="Z") +dist_speed <- do_fit_velocity(data$speed, suffix="vtb") +dist_direction <- do_fit_direction(data$direction[data$direction<=0], suffix="vtb") +model_acf_x <- do_fit_acf(data$wind_velocity_x, axis="X", suffix="vtb") +model_acf_y <- do_fit_acf(data$wind_velocity_y, axis="Y", suffix="vtb") +model_acf_z <- do_fit_acf(data$wind_velocity_z, axis="Z", suffix="vtb") banner("Anemometer") data_real <- select_samples(time_start, time_delta) -dist_speed_real <- do_fit_velocity(data_real$speed) -dist_direction_real <- do_fit_direction(data_real$direction) -model_acf_x <- do_fit_acf(data_real$x[data_real$x<=0], axis="X") -model_acf_y <- do_fit_acf(data_real$y[data_real$y<=0], axis="Y") -model_acf_z <- do_fit_acf(data_real$z[data_real$z<=0], axis="Z") +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") +model_acf_y_real <- do_fit_acf(data_real$y[data_real$y<=0], axis="Y", suffix="anem") +model_acf_z_real <- do_fit_acf(data_real$z[data_real$z<=0], axis="Z", suffix="anem") -#print(data) +# 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 +} + +banner("Comparison") +plot(dist_speed$x, normalise(dist_speed$actual_density), xlab='Wind speed', ylab='PDF') +lines(dist_speed$x, normalise(dist_speed$estimated_density), col='red', lwd=2) +points(dist_speed_real$x, normalise(dist_speed_real$actual_density), col='blue') +lines(dist_speed_real$x, normalise(dist_speed_real$estimated_density), col='blue', lwd=2) + +plot(dist_direction$x, normalise(dist_direction$actual_density), xlab='Wind direction', ylab='PDF') +lines(dist_direction$x, normalise(dist_direction$estimated_density), col='red', lwd=2) +points(dist_direction_real$x, normalise(dist_direction_real$actual_density), col='blue') +lines(dist_direction_real$x, normalise(dist_direction_real$estimated_density), col='blue', lwd=2) + +plot(model_acf_x$lag, model_acf_x$actual, col='red', xlab='Lag', ylab="ACF X") +lines(model_acf_x$lag, model_acf_x$estimated, col='red', lwd=2) +points(model_acf_x_real$lag, model_acf_x_real$actual, col='blue') +lines(model_acf_x_real$lag, model_acf_x_real$estimated, col='blue', lwd=2) + +plot(model_acf_y$lag, model_acf_y$actual, col='red', xlab='Lag', ylab="ACF Y") +lines(model_acf_y$lag, model_acf_y$estimated, col='red', lwd=2) +points(model_acf_y_real$lag, model_acf_y_real$actual, col='blue') +lines(model_acf_y_real$lag, model_acf_y_real$estimated, col='blue', lwd=2) + +plot(model_acf_z$lag, model_acf_z$actual, col='red', xlab='Lag', ylab="ACF Z") +lines(model_acf_z$lag, model_acf_z$estimated, col='red', lwd=2) +points(model_acf_z_real$lag, model_acf_z_real$actual, col='blue') +lines(model_acf_z_real$lag, model_acf_z_real$estimated, col='blue', lwd=2) diff --git a/gnuplot/vtestbed.gnuplot b/gnuplot/vtestbed.gnuplot @@ -0,0 +1,66 @@ +set terminal svg size 450,450/2 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 title offset 0,-1 +set tmargin 1 +set lmargin 6 +set rmargin 1 +unset key + +set output 'build/gnuplot/vtestbed.svg' + +set multiplot layout 2,3 + +n = 1 +set xrange [0:40] +set yrange [0:*] +set xlabel 'Lag' + +set title 'ACF X' +set xtics 0,10 +plot \ +'build/vtestbed/anem/acf-X' using 1:2 every n with points pt 6 lc '#404040' notitle,\ +'' using 1:3 with lines lw 2 lc '#404040' title 'Real',\ +'build/vtestbed/vtb/acf-X' using 1:2 every n with points pt 6 lc '#c04040' notitle,\ +'' using 1:3 with lines lw 2 lc '#c04040' title 'Simulated' + +set title 'ACF Y' +plot \ +'build/vtestbed/anem/acf-Y' using 1:2 every n with points pt 6 lc '#404040' notitle,\ +'' using 1:3 with lines lw 2 lc '#404040' title 'Real',\ +'build/vtestbed/vtb/acf-Y' using 1:2 every n with points pt 6 lc '#c04040' notitle,\ +'' using 1:3 with lines lw 2 lc '#c04040' title 'Simulated' + +set title 'ACF Y' +plot \ +'build/vtestbed/anem/acf-Z' using 1:2 every n with points pt 6 lc '#404040' notitle,\ +'' using 1:3 with lines lw 2 lc '#404040' title 'Real',\ +'build/vtestbed/vtb/acf-Z' using 1:2 every n with points pt 6 lc '#c04040' notitle,\ +'' using 1:3 with lines lw 2 lc '#c04040' title 'Simulated' + +n = 5 + +set xrange [0:*] +set xtics 0,2 +set xlabel 'Velocity, m/s' +set title 'Wind speed' +plot \ +'build/vtestbed/anem/velocity' using 1:2 every n with points pt 6 lc '#404040' notitle,\ +'' using 1:3 with lines lw 2 lc '#404040' title 'Real',\ +'build/vtestbed/vtb/velocity' using 1:2 every n with points pt 6 lc '#c04040' notitle,\ +'' using 1:3 with lines lw 2 lc '#c04040' title 'Simulated' + +n = 20 + +set xrange [*:0] +set xtics -4,1 +set key top right Left reverse width -20 +set xlabel 'Direction, degrees' +set title 'Wind direction' +plot \ +'build/vtestbed/anem/direction' using 1:2 every n with points pt 6 lc '#404040' notitle,\ +'' using 1:3 with lines lw 2 lc '#404040' title 'Real',\ +'build/vtestbed/vtb/direction' using 1:2 every n with points pt 6 lc '#c04040' notitle,\ +'' using 1:3 with lines lw 2 lc '#c04040' title 'Simulated' diff --git a/main.bib b/main.bib @@ -160,3 +160,16 @@ year = {1988}, institution = {Sandia National Labs., Albuquerque, NM (USA)} } + +@InProceedings{ bogdanov2020vtestbed, + author = {Bogdanov, Alexander and Degtyarev, Alexander and Gankevich, + Ivan and Khramushin, Vasily and Korkhov, Vladimir}, + title = {Virtual Testbed: Concept and Applications}, + booktitle = {Computational Science and Its Applications -- ICCSA 2020}, + year = {2020}, + publisher = {Springer International Publishing}, + address = {Cham}, + pages = {3--17}, + isbn = {978-3-030-58817-5}, + doi = {10.1007/978-3-030-58817-5\_1} +} diff --git a/main.tex b/main.tex @@ -23,7 +23,7 @@ \lstset{language=R,literate={<-}{{$\gets$}}1} -\titlerunning{TODO} +\titlerunning{Wind simulation using high-frequency velocity component measurements} \authorrunning{A.\,Gavrikov, I.\,Gankevich} \institute{Saint Petersburg State University\\ @@ -43,7 +43,8 @@ \and strain gauge \and anemometer \and three-dimensional ACF - \and wind velocity PDF. + \and wind velocity PDF + \and autoregressive model. } \end{abstract} @@ -369,7 +370,7 @@ sampleToSpeed <- function(x, c1, c2) { Over a period of one month we collected 3.1M samples and filtered out 12\% of them as having too large unnatural values. We attributed these values to measurement errors as they spread uniformly across all the time span and are -surrounded by the values of regular magnitude. TODO picture of a large error. +surrounded by the values of regular magnitude. From the remaning data we choose days with wind speeds above the average as reported by EMERCOM of Russia\footnote{https://en.mchs.gov.ru/} for Saint @@ -454,10 +455,6 @@ similar shapes, for \(x\) and \(y\) axes the distribution for incident flow is taller than the distribution for turbulent flow. The best-fit and worst-fit distributions for each axis are presented in 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 \includegraphics{build/gnuplot/velocity-dist.eps} @@ -521,8 +518,6 @@ correspondence between the measurements of the two anemometers anemometers.\label{fig-hold-peak}} \end{figure} -%TODO asymmetry - Finally, we computed autocovariance for each axis as \begin{equation} K(\tau) = \Expectation\left[ \left(X_t-\bar{X}\right)\left(X_{t-\tau}-\bar{X}\right) \right] @@ -551,7 +546,40 @@ are presented in figure~\ref{fig-acf}. \subsection{Wind simulation using measured ACFs} -TODO Anton +We simulated three-dimensional wind velocity using autoregressive model +implemented in Virtual Testbed~\cite{bogdanov2020vtestbed} and using the data +obtained with three-axis anemometer on March 28, 2021, 01:00-03:00 UTC. We +approximated four-dimensional autocovariance function using~\eqref{eq-acf} by +setting the corresponding parameters from one-dimensional autocovariance +functions estimated from the data obtained with anemometer, all other +parameters were set close to nought. Non-nought parameters are listed in +table~\ref{tab-arma-parameters}. We found that velocity and direction +distributions and ACFs of each axis of simulated wind and real wind are close +to each other (see~figure~\ref{fig-vtestbed}). + +\begin{table} + \centering + \caption{Input parameters for AR model that were used to simulate wind velocity. + \label{tab-arma-parameters}} + \begin{tabular}{lllll} + \toprule + Axis & ACF \(a\) & \phantom{0}ACF \(b\) & \phantom{0}ACF \(c\) & \phantom{0}Mean velocity \\ + \midrule + \(x\) & 1.793 & \phantom{0}0.0214 & \phantom{0}0.2603 & \phantom{0}-2.439 \\ + \(y\) & 1.423 & \phantom{0}0.01429 & \phantom{0}0.2852 & \phantom{0}-2.158 \\ + \(z\) & 0.9075 & \phantom{0}0.06322 & \phantom{0}0.3349 & \phantom{0}-1.367 \\ + \bottomrule + \end{tabular} +\end{table} + +\begin{figure} + \centering + \includegraphics{build/gnuplot/vtestbed.eps} + \caption{Comparison of simulated wind velocity and direction distributions + and per-axis ACFs + between the data from the anemometer and the data from Virtual Testbed + (2021, March 28, 01:00--03:00, UTC).\label{fig-vtestbed}} +\end{figure} \section{Discussion}