iccsa-21-wind

Wind Simulation Using High-Frequency Velocity Component Measurements
git clone https://git.igankevich.com/iccsa-21-wind.git
Log | Files | Refs

commit 4b9052e50f2e317d42da28fe2e65e2a9f7e3fac5
parent d55b7cb3c45107ebe011f0a6bfedee7da046a92d
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date:   Sat,  1 Jan 2022 22:10:07 +0300

Merge remote-tracking branch 'refs/remotes/origin/master'

Diffstat:
.gitignore | 1+
Makefile | 3+--
copyright.txt | 2++
main.tex | 65++++++++++++++++++++++++++++++++++++++++-------------------------
photos/anemometer.jpg | 0
slides.tex | 261+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
6 files changed, 301 insertions(+), 31 deletions(-)

diff --git a/.gitignore b/.gitignore @@ -1,5 +1,6 @@ /build Rplots.pdf +/CoA_Medium.eps # Created by https://www.toptal.com/developers/gitignore/api/vim,linux,r # Edit at https://www.toptal.com/developers/gitignore?templates=vim,linux,r diff --git a/Makefile b/Makefile @@ -72,8 +72,7 @@ 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 +build/main.zip: main.tex main.bib llncs.cls splncs04.bst build/inkscape/anemometer.eps build/daily-stats.eps build/gnuplot/velocity-dist.eps build/gnuplot/velocity-xy-dist-max.eps build/gnuplot/velocity-xy-dist-min.eps build/gnuplot/direction-dist.eps build/gnuplot/hold-peak.eps build/gnuplot/acf.eps build/gnuplot/turbulence.eps build/gnuplot/vtestbed.eps @mkdir -p build zip --filesync build/main.zip $^ diff --git a/copyright.txt b/copyright.txt @@ -4,3 +4,5 @@ Anton Gavrikov, Ivan Gankevich Wind simulation in the context of ship motions is used to estimate the effect of the wind on large containerships, sailboats and yachts. Wind models are typically based on a sum of harmonics with random phases and different amplitudes. In this paper we propose to use autoregressive model to simulate the wind. This model is based on autocovariance function that can be estimated from the real-world data collected by anemometers. We have found none of the data that meets our resolution requirements, and decided to produce the dataset ourselves using three-axis anemometer. We built our own anemometer based on load cells, collected the data with the required resolution, verified the data using well-established statistical distributions, estimated autocovariance functions from the data and simulated the wind using autoregressive model. We have found that the load cell anemometer is capable of recording wind speed for statistical studies, but autoregressive model needs further calibration to reproduce the wind with the same statistical properties. Keywords: load cell, strain gauge, anemometer, three-dimensional ACF, wind velocity PDF, autoregressive model, turbulence. + +Ivan Gankevich, Saint Petersburg State University, 13B Universitetskaya Emb., St Petersburg 199034, Russia, i.gankevich@spbu.ru diff --git a/main.tex b/main.tex @@ -8,6 +8,7 @@ \usepackage{listings} \usepackage{url} \usepackage{textcomp} +\usepackage{xcolor} \DeclareMathOperator{\Mode}{Mo} \DeclareMathOperator{\Expectation}{E} @@ -21,7 +22,19 @@ Ivan Gankevich\textsuperscript{*}\orcidID{0000-0001-7067-6928} } -\lstset{language=R,literate={<-}{{$\gets$}}1} +\lstdefinelanguage{rrr}{ + basicstyle={\rmfamily}, + morekeywords={function}, + keywordstyle={\bfseries}, + showstringspaces=false, + texcl=true, + escapeinside={LATEX}{END}, + otherkeywords={function} +} +\lstset{ + language=rrr, + literate={<-}{{$\gets$}}1, +} \titlerunning{Wind simulation using high-frequency velocity component measurements} \authorrunning{A.\,Gavrikov, I.\,Gankevich} @@ -97,7 +110,7 @@ coefficients are different~\cite{veers1984modeling}. %& c_1 = 4 && c_2 = 70 && \text{for }y,\\ %& c_1 = 0.5 && c_2 = 8 && \text{for }z. %\end{align*} -The spectrum describes wind velocity vector in the plane that is perpendicular +The spectrum describes wind velocity vector in the plane that is perpendicular to the mean wind direction vector and travels in the same direction with mean wind speed. Time series is generated as Fourier series, coefficients of which are determined from the spectrum, and phases are random @@ -136,7 +149,7 @@ a decaying exponent: \end{equation*} where \(\sigma^2\) is process variance (area under the spectrum). -Unfortunately, this autocovariance function is one-dimensional and there is +Unfortunately, this autocovariance function is one-dimensional and there is no easy way of obtaining three-dimensional analogue from the spectra. In order to solve this problem we looked into datasets of wind speed measurements available for the research. However, most of them contain either one or two @@ -186,7 +199,7 @@ three axes at one point in space multiple times per second. To measure wind velocity we used resistive foil strain gauges mounted on the arms aligned perpendicular to the axes directions. Inside each arm we placed aluminium load cell with two strain gauges: one on the bending side and another -one on the lateral side. Load cells use Wheatstone bridge to measure the +one on the lateral side. Load cells use Wheatstone half-bridge to measure the resistance of the gauges and are connected to the circuit that measures the resistance and transmits it to the microcontroller in digital format. Microcontroller then records the value for each load cell and transmits all of @@ -210,8 +223,9 @@ figure~\ref{fig-anemometer}. \end{tabular} \end{minipage} %\includegraphics{build/gnuplot/anemometer.eps} - \caption{Three-axis anemometer. Arms are inserted into the Housing - and fixed using bolts that go through the circular holes.\label{fig-anemometer}} + \caption{Three-axis anemometer. Arms are inserted into the housing + and fixed using bolts that go through the circular holes. Red, green, blue colours +denote load cells for \(x\), \(y\), \(z\) axes respectively.\label{fig-anemometer}} \end{figure} Each load cell faces the direction that is perpendicular to the directions of @@ -249,7 +263,7 @@ For all load cells we have a system of three equations Hence \(\upsilon_{x}=\alpha_{x}\sqrt{F_{x}}\), \(\upsilon_{y}=\alpha_{y}\sqrt{F_{y}}\), \(\upsilon_{z}=\alpha_{z}\sqrt{F_{z}}\) where \(\alpha_{x,y,z}\) are constants of -proportionality. Therefore, to obtain wind velocity we \emph{take square root} of the +proportionality. Therefore, to obtain wind velocity we take \emph{square root} of the value measured by the load cell and multiply it by some coefficient determined empirically during anemometer calibration. @@ -267,7 +281,7 @@ scalar wind velocity is a length of wind velocity vector \(\upsilon=\sqrt{\smash[b]{\upsilon_x^2+\upsilon_y^2+\upsilon_z^2}}\). However, projection of wind velocity vector on \(x\), \(y\) or \(z\) axis may be negative. Our solution to this problem is to use two Weibull distributions: one for positive and -one for negative projection --- but with different parameters. +one for negative projection~--- but with different parameters. \begin{equation} f\left(\upsilon_x;b_1,c_1,b_2,c_2\right) = \begin{cases} b_1 c_1 \left(b_1 \left|\upsilon_x\right|\right)^{c_1-1} \exp\left(-\left(b_1\left|\upsilon_x\right|\right)^{c_1}\right). @@ -357,8 +371,8 @@ placed the fan behind and in front of each arm of our anemometer and measured values that the corresponding load cell reported for each side of the arm with and without the fan. Then we were able to calculate two coefficients for each axis: one for negative and another one for positive wind velocities. The -coefficient equals the raw sensor value that is equivalent to the wind speed of 1 -m/s (table~\ref{tab-coefficients}). +coefficient equals the raw sensor value that is equivalent to the wind speed of 1~m/s +(table~\ref{tab-coefficients}). \begin{table} \begin{minipage}[t]{0.35\textwidth} @@ -376,6 +390,7 @@ m/s (table~\ref{tab-coefficients}). \(C_1\) is for negative values and \(C_2\) is for positive values. \label{tab-coefficients}} \end{minipage} + ~ \begin{minipage}[t]{0.55\textwidth} \centering \begin{tabular}{ll} @@ -398,20 +413,20 @@ removed this linear trend from the measured values using linear regression. The code in R~\cite{r-language} that transforms raw sensor values into wind speed is presented in listing~\ref{lst-sample-to-speed}. -\begin{lstlisting}[label={lst-sample-to-speed},caption={The code that transforms raw load cell raw sensor values into wind speed projections to the corresponding axis.}] +\begin{lstlisting}[label={lst-sample-to-speed},caption={The code that transforms raw load cell sensor values into wind speed projections to the corresponding axis.}] sampleToSpeed <- function(x, c1, c2) { t <- c(1:length(x)) reg <- lm(x~t) - x <- x - reg$fitted.values # remove linear trend - x <- sign(x)*sqrt(abs(x)) # convert from force to velocity - x[x<0] = x[x<0] / c1 # scale sensor values to wind speed - x[x>0] = x[x>0] / c2 # using calibration coefficients + x <- x - reg$fitted.values LATEX{\color{black!70}\textit{\# remove linear trend}}END + x <- sign(x)*sqrt(abs(x)) LATEX{\color{black!70}\textit{\# convert from force to velocity}}END + x[x<0] = x[x<0] / c1 LATEX{\color{black!70}\textit{\# scale sensor values to wind speed}}END + x[x>0] = x[x>0] / c2 LATEX{\color{black!70}\textit{\# using calibration coefficients}}END x } \end{lstlisting} 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 +them because they had 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. %From the remaining data we choose days with wind speeds above the average as @@ -427,9 +442,9 @@ figure~\ref{fig-intervals}, dataset properties are presented in table~\ref{tab-d Unique feature of three-axis anemometer is that it measures both velocity of incident air flow towards the arms and the turbulent flow that forms behind the -arms. Turbulent flow velocity distribution mode is often smaller than the mode +arms. Turbulent flow velocity distribution peak is often smaller than the peak of incident flow. To exclude it from the measurements one can choose the side -with the largest mode (either positive or negative part +with the largest peak (either positive or negative part of~\eqref{eq-velocity-distribution}), but in this paper we left them as is for the purpose of the research. @@ -488,9 +503,9 @@ distributions for each axis are presented in figure~\ref{fig-velocity-distributi Wind direction was approximated by von Mises distribution using least-squares fitting. Following the common practice~\cite{feng2015sectors} we divided direction axis -into sectors: from -180\textdegree{} to -90\textdegree{}, from from -90\textdegree{} to +into sectors: from -180\textdegree{} to -90\textdegree{}, from -90\textdegree{} to 0\textdegree{}, from 0 to 90\textdegree{} and from 90\textdegree{} to -180\textdegree{} --- and fitted each sector independently. We chose four sector +180\textdegree{}~--- and fitted each sector independently. We chose four sector to have one sector for each side of the anemometer. The best-fit and worst-fit distributions for each sector are presented in figure~\ref{fig-direction-distributions}. @@ -521,7 +536,7 @@ are presented in figure~\ref{fig-direction-distributions}. In order to verify that our anemometer produces correct time series we performed synchronous measurements with commercial anemometer HP-866A. Two anemometers were placed in the open field near the shore. The commercial -anemometer directed towards the mean wind direction as it measures the speed +anemometer was directed towards the mean wind direction as it measures the speed only. The data from the two anemometers was recorded synchronously. From the data we selected intervals with wind gusts and compared the measurements of the two anemometers. To compare them we scaled load cell anemometer measurements to @@ -559,7 +574,7 @@ are presented in figure~\ref{fig-acf}. data intervals with the smallest error (the second row). Red line shows estimated ACF of positive wind speed projections, blue line shows estimated ACF of negative wind speed projections and circles denote observed - ACF of wind speed projections. + ACF of wind speed projections. Note: graphs use mirrored axes. \label{fig-acf}} \end{figure} @@ -600,7 +615,7 @@ 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 similar in shape, but are too far away from each other (see~figure~\ref{fig-vtestbed}). We consider these results preliminary and -will investigate further in future work. +will investigate them further in future work. \begin{table} \centering @@ -673,8 +688,8 @@ largest NRMSE of 6.7\% and wind directions fit into von Mises distribution with the largest NRMSE of 11\%. We estimated autocovariance functions for wind speed for each axis of the anemometer and used this approximations to simulate wind flow in Virtual Testbed. The parameters of these functions allow to control -both wind speed and mean direction. The future work is to construct new -anemometer that is able to measure spatial autocovariance using the proposed +both wind speed and mean direction. The future work is to construct an array of +anemometers that is able to measure spatial autocovariance using the proposed anemometer as the base. \subsubsection*{Acknowledgements.} diff --git a/photos/anemometer.jpg b/photos/anemometer.jpg Binary files differ. diff --git a/slides.tex b/slides.tex @@ -1,7 +1,7 @@ \documentclass[aspectratio=169,xcolor=table]{beamer} \usepackage{polyglossia} \setdefaultlanguage{english} -\usetheme{SaintPetersburg} +\usetheme[numbers]{SaintPetersburg} \usepackage{textcomp} \usepackage{booktabs} @@ -24,22 +24,275 @@ pdfmetalang={en}, pdflicenseurl={http://creativecommons.org/licenses/by-sa/4.0/}, pdfcopyright={Copyright \textcopyright{} 2021 Anton Gavrikov\xmpcomma{} Ivan Gankevich\xmpcomma{}}, - pdfsubject={TODO}, + pdfsubject={Wind simulation using high-frequency velocity component measurements}, } -\title{TODO} +\title{Wind simulation using high-frequency velocity component measurements} \author{Anton Gavrikov \and Ivan Gankevich} \institute{Saint Petersburg State University} -\date{July 2021} +\date{September 2021} + +\lstdefinelanguage{rrr}{ + basicstyle={\rmfamily\footnotesize}, + morekeywords={function}, + keywordstyle={\bfseries}, + showstringspaces=false, + texcl=true, + escapeinside={LATEX}{END}, + otherkeywords={function} +} +\lstset{ + language=rrr, + literate={<-}{{$\gets$}}1, +} \begin{document} \frame{\maketitle} \begin{frame}{Motivation} + \begin{itemize} + \item Datasets with all three wind velocity components and one-second resolution are difficult to find. + \item Three-axis ultrasonic anemometers are expensive. + \item We needed three-dimensional autocovariance function to simulate wind with autoregressive + model. + \end{itemize} + \vfill + Solution: make our own anemometer based on load cells. +\end{frame} + +\begin{frame}{Three-axis anemometer based on load cells} + \centering + \begin{columns} + \begin{column}{0.45\textwidth} + \includegraphics{build/inkscape/anemometer.eps} + \end{column} + \begin{column}{0.45\textwidth} + \includegraphics[width=0.85\textwidth]{photos/anemometer.jpg} + \vskip0.25\baselineskip + \begin{tabular}{ll} + Load cell capacity & 1 kg \\ + Load cell amplifier & HX711 \\ + Microcontroller & ATmega328P \\ + \end{tabular} + \end{column} + \end{columns} +\end{frame} + +\begin{frame}{Calculating wind velocity} + Bernoulli's equation: + \begin{equation*} + \rho\frac{\upsilon^2}{2} + \rho gz = p_0-p + \end{equation*} + Pressure force: + \begin{equation*} + \vec{F} = p S \vec{n}. + \end{equation*} + Wind velocity: + \begin{equation*} + \begin{cases} + \upsilon_x^2 \propto F_x \\ + \upsilon_y^2 \propto F_y \\ + \upsilon_z^2 \propto F_z \\ + \end{cases} \Rightarrow + \begin{cases} + \upsilon_{x}=\alpha_{x}\sqrt{F_{x}} \\ + \upsilon_{y}=\alpha_{y}\sqrt{F_{y}} \\ + \upsilon_{z}=\alpha_{z}\sqrt{F_{z}} \\ + \end{cases} + \end{equation*} +\end{frame} + +\begin{frame}{Per-axis probability distribution function for wind velocity} + Weibull distribution (scalar velocity)\footnote{C. Justus, W. Hargraves, A. Yalcin + \textit{Nationwide Assessment of Potential Output from Wind-Powered Generators}, 1976.}: + \begin{equation*} + f\left(\upsilon;b,c\right) = b c \left(b \upsilon\right)^{c-1} + \exp\left(-\left(b\upsilon\right)^{c}\right). + \end{equation*} + Weibull distribution (per-axis scalar velocity): + \begin{equation*} + f\left(\upsilon_x;b_1,c_1,b_2,c_2\right) = \begin{cases} + b_1 c_1 \left(b_1 \left|\upsilon_x\right|\right)^{c_1-1} \exp\left(-\left(b_1\left|\upsilon_x\right|\right)^{c_1}\right). + \qquad \text{if } \upsilon_x < 0 \\ + b_2 c_2 \left(b_2 \left|\upsilon_x\right|\right)^{c_2-1} \exp\left(-\left(b_2\left|\upsilon_x\right|\right)^{c_2}\right). + \qquad \text{if } \upsilon_x \geq 0 \\ + \end{cases} + \end{equation*} + Parameters: + + \vspace{0.5\baselineskip} + \begin{tabular}{ll} + \(b\) & scale parameter \\ + \(c\) & shape parameter \\ + \end{tabular} +\end{frame} + +\begin{frame}{Three-dimensional ACF of wind velocity} + One-dimensional autocovariance function\footnote{G. Box, G. Jenkins \textit{Time series analysis: forecasting and control}, 1976.}: + \begin{equation*} + K\left(t\right) = a_3 \exp\left(-\left(b_3 t\right)^{c_3} \right). + \end{equation*} + Three-dimensional autocovariance function: + \begin{equation*} + K\left(t,x,y,z\right) = a \exp\left( + -\left(b_t t\right)^{c_t} + -\left(b_x x\right)^{c_x} + -\left(b_y y\right)^{c_y} + -\left(b_z z\right)^{c_z} + \right). + \end{equation*} +\end{frame} + +\begin{frame}[fragile]{Data collection and preprocessing} + \begin{columns}[T] + \begin{column}{0.30\textwidth} + Calibration coefficients:\strut{} + \begin{tabular}{lll} + \toprule + Axis & \(C_1\) & \(C_2\) \\ + \midrule + X & 11.19 & 12.31 \\ + Y & 11.46 & 11.25 \\ + Z & 13.55 & 13.90 \\ + \bottomrule + \end{tabular} + \end{column} + \begin{column}{0.65\textwidth} + Dataset properties:\strut{} + \begin{tabular}{ll} + \toprule + Time span & 36 days \\ + Size & 122 Mb \\ + No. of samples & 3\,157\,234 \\ + No. of samples after filtering & 2\,775\,387 \\ + Resolution & 1 sample per second \\ + \bottomrule + \end{tabular} + \end{column} + \end{columns} + \vskip0.5\baselineskip + \(C_1\) --- negative values, \(C_2\) --- positive values. + \vskip0.25\baselineskip +\begin{lstlisting} +sampleToSpeed <- function(x, c1, c2) { + t <- c(1:length(x)) + reg <- lm(x~t) + x <- x - reg$fitted.values LATEX{\color{black!70}\textit{\# remove linear trend}}END + x <- sign(x)*sqrt(abs(x)) LATEX{\color{black!70}\textit{\# convert from force to velocity}}END + x[x<0] = x[x<0] / c1 LATEX{\color{black!70}\textit{\# scale sensor values to wind speed}}END + x[x>0] = x[x>0] / c2 LATEX{\color{black!70}\textit{\# using calibration coefficients}}END + x } +\end{lstlisting} +\end{frame} + +\begin{frame}{Verification against EMERCOM data} + \begin{columns} + \begin{column}{0.70\textwidth} + \includegraphics[height=0.9\textheight]{build/daily-stats.eps} + \end{column} + \begin{column}{0.25\textwidth} + \footnotesize + \(\upsilon_x\), \(\upsilon_y\), \(\upsilon_z\), \(\upsilon\) --- mean velocities. + + \vskip\baselineskip + Yellow regions --- wind velocity above average as reported by EMERCOM of Russia + \href{https://en.mchs.gov.ru/}{https://en.mchs.gov.ru/}. + \end{column} + \end{columns} +\end{frame} + +\begin{frame}{Verification of wind velocity against Weibull distribution} + \begin{columns}[T] + \begin{column}{0.80\textwidth} + \includegraphics{build/gnuplot/velocity-dist.eps} + \end{column} + \begin{column}{0.15\textwidth} + \vskip2\baselineskip + largest error + \vskip6\baselineskip + smallest error + \end{column} + \end{columns} +\end{frame} + +\begin{frame}{Verification of wind direction against von Mises distribution} + \centering + \includegraphics{build/gnuplot/direction-dist.eps} + + \phantom{mm}largest error\phantom{mmmmmmmmmm}smallest error + \vfill + \begin{flushleft} + \footnotesize{} + J. Carta, C. Bueno, P. Ramírez + \textit{Statistical modelling of directional wind speeds using mixtures of + von Mises distributions: Case study}, 2008. + \end{flushleft} +\end{frame} + +\begin{frame}{Verification of wind velocity against commercial anemometer} + \centering + \includegraphics{build/gnuplot/hold-peak.eps} + + \phantom{mm}first site\phantom{mmmmmmmmmm}second site + \footnotesize + %\begin{equation} + % \text{NRMSE} = \frac{\sqrt{\Expectation\left[\left(X_\text{observed}-X_\text{estimated}\right)^2\right]}}{X_\text{max}-X_\text{min}}. + %\end{equation} +\end{frame} + +\begin{frame}{Verification of ACF against approximation} + \begin{columns}[T] + \begin{column}{0.80\textwidth} + \includegraphics{build/gnuplot/acf.eps} + \end{column} + \begin{column}{0.15\textwidth} + \vskip2\baselineskip + largest error + \vskip6\baselineskip + smallest error + \end{column} + \end{columns} +\end{frame} + +\begin{frame}{Turbulence coefficient} + The ratio of absolute mean speed of turbulent flow to the absolute mean + speed of incident flow for each axis. + \begin{center} + \includegraphics{build/gnuplot/turbulence.eps} + \end{center} +\end{frame} + +\begin{frame}{Wind simulation} + \includegraphics{build/gnuplot/vtestbed.eps} + \begin{tikzpicture}[overlay] + \node at (-2.50,1) {\footnotesize{}Model parameters:}; + \node at (-0.25,-0.20) { + \footnotesize + \begin{tabular}{lllll} + \toprule + & ACF \(a\) & ACF \(b\) & ACF \(c\) & Mean \(\upsilon\), m/s \\ + \midrule + \(x\) & 1.793 & 0.0214 & 0.2603 & -2.439 \\ + \(y\) & 1.423 & 0.01429 & 0.2852 & -2.158 \\ + \(z\) & 0.9075 & 0.06322 & 0.3349 & -1.367 \\ + \bottomrule + \end{tabular} + }; + \end{tikzpicture} \end{frame} \begin{frame}{Conclusion and future work} + \begin{itemize} + \item Per-axis wind speeds fit into Weibull distribution with max. NRMSE of 6.7\%. + \item Wind directions fit into von Mises distribution with max. NRMSE of 11\%. + \item Three-axis anemometer is useful for statistical studies, + but is unable to measure \textit{immediate} wind speed and direction. + \item Simulated with autoregressive model and real wind velocity are + similar in shape, but too far away from each other. + \end{itemize} + \vfill + Future work: construct an array of anemometers to measure spatial autocovariance. \end{frame} \begin{frame}