1 Wczytanie danych i wykresy

Zaczniemy od wczytania pliku Excel dane_w2.xlsx (być może trzeba ściągnąć ten plik ze strony przedmiotu jeszcze raz lub podać odpowiednią lokalizację pliku na dysku. Można wykorzystać też kreator importu):

library(readxl)
dane_w2 <- read_excel("dane_w2.xlsx")
head(dane_w2)
## # A tibble: 6 × 49
##   `Date & Time Stamp` CH1Avg CH1SD CH1Max CH1Min CH2Avg CH2SD CH2Max CH2Min
##   <dttm>               <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>
## 1 2010-08-01 00:00:00    4.6   0.3    5.3    3.8      0     0      0      0
## 2 2010-08-01 00:10:00    4.7   0.2    5.3    4.1      0     0      0      0
## 3 2010-08-01 00:20:00    5.1   0.3    5.7    4.5      0     0      0      0
## 4 2010-08-01 00:30:00    5.4   0.3    6      4.9      0     0      0      0
## 5 2010-08-01 00:40:00    5     0.3    5.7    4.1      0     0      0      0
## 6 2010-08-01 00:50:00    4.1   0.2    4.5    3.4      0     0      0      0
## # ℹ 40 more variables: CH3Avg <dbl>, CH3SD <dbl>, CH3Max <dbl>, CH3Min <dbl>,
## #   CH4Avg <dbl>, CH4SD <dbl>, CH4Max <dbl>, CH4Min <dbl>, CH5Avg <dbl>,
## #   CH5SD <dbl>, CH5Max <dbl>, CH5Min <dbl>, CH6Avg <dbl>, CH6SD <dbl>,
## #   CH6Max <dbl>, CH6Min <dbl>, CH7Avg <dbl>, CH7SD <dbl>, CH7Max <dbl>,
## #   CH7Min <dbl>, CH8Avg <dbl>, CH8SD <dbl>, CH8Max <dbl>, CH8Min <dbl>,
## #   CH9Avg <dbl>, CH9SD <dbl>, CH9Max <dbl>, CH9Min <dbl>, CH10Avg <dbl>,
## #   CH10SD <dbl>, CH10Max <dbl>, CH10Min <dbl>, CH11Avg <dbl>, CH11SD <dbl>, …

Istotne składniki zestawu pomiarowego:

  1. Anemometry (prędkość wiatru):
    CH1, 81,00m; (backup - zapasowy)
    CH4, 80,60m;
    CH5, 60,90m;
    CH6, 40,90m

  2. Vane’y (kierunek wiatru):
    CH7, 78,30m;
    CH8, 38,60m

  3. Inne:
    CH9 wilgotność
    CH10 ciśnienie
    CH11 temperatura
    CH12 napięcie w rejestratorze

Końcówki w nazwach kolumn odnoszą się do:

1.1 Zadanie - wykres prędkości wiatru

Narysować wykres (ggplot2) zależności prędkości wiatru dla wszystkich anemometrów (zmienne CH1Avg, CH4Avg, CH5Avg, CH6Avg) od czasu na jednym wykresie. Nadać wykresowi tytuł, opisy osi (w języku polskim).

Wykres powinien przypominać:

1.2 Zadanie - ggplotly

Powyższy wykres ze względu na dużą ilość danych może być mało czytelny. Przy użyciu pakietu plotly można zrobić jego interaktywną wersję - link do przykładu

Wskazówka: należy wykorzystać rysunek z poprzedniego zadania i przekazać go do odpowiedniej funkcji pakietu plotly.

Uwaga: możliwe, że w laboratoriach wydziału jest stara wersja pakietu plotly i przykład z linka może nie działać. Należy zainstalować ten pakiet ponownie (lub zaktualizować). Z przykładu potrzebne nam są tylko biblioteki ggplot2 oraz plotly.

Przykładowy wynik:

1.3 Zadanie - interakcja z wykresem

Na stworzonym wykresie wypróbować:

  • powiększenie wybranego fragmentu wykresu

  • przesuwanie wykresu

  • włączenie/wyłączenie serii danych poprzez kliknięcie ich na legendzie.

1.4 Interpretacja wykresów

Wykresy prędkości wiatru w czasie pozwalają na wykrycie nieprawidłowości w działaniu urządzeń pomiarowych:

  • Anemometry CH1 i CH4 zostały zamontowane na tej samej wysokości, w związku z czym ich wskazania powinny być zbliżone (niemal identyczne)

  • Wraz ze wzrostem wysokości (na ogół) wzrasta prędkość wiatru, zatem jeżeli wysokości zamontowanych urządzeń to:
    CH1, 81.00m; (backup - zapasowy)
    CH4, 80.60m;
    CH5, 60.90m;
    CH6, 40.90m
    to wykresy powinny układać się “monotonicznie”. Co prawda, mogą zdarzać się inwersje np. przy niskich wartościach, jednak może być to wskazówka o awarii urządzenia

  • W miesiącach zimowych wykres jak poniżej pozwala na szybkie zidentyfikowanie sytuacji, gdy urządzenia pomiarowe zamarzły - płaska linia na poziomie wartości 0 - tutaj widoczne problemy w końcówce listopada:

2 Profil pionowy wiatru

2.1 Wzór/wykładnik Hellmana-Suttona.

Prędkość (średnia) wiatru nie zależy od wysokości w sposób liniowy, jednak podlega wzorowi Hellmana-Suttona:

\[ V(h) = V_1 \cdot \left( \frac{h}{h_1} \right)^w \quad (1) \]

gdzie:
\(h_1\) - bazowa wysokość - znamy wysokości zamontowanych urządzeń, wybieramy jedno z nich
\(V_1\) - średnia prędkość zmierzona na bazowej wysokości - jesteśmy w stanie ją wyliczyć z danych pomiarowych
\(w\) - “nieznany”, ale możliwy do wyznaczenia wykładnik. Jest on charakterystyczny dla danego miejsca pomiaru.

Znając powyższe 3 wartości (po wyliczeniu są to stałe) dostajemy funkcję jednej zmiennej opisującą zależność prędkości wiatru od wysokości.

2.2 Zadanie - wzór Hellmana-Suttona

Przyjmijmy, że dane \(h_1\) i \(V_1\) dotyczą anemometru CH6 oraz znając wartości \(h_2\) (z danych o instalacji) i \(V_2\) (trzeba wyliczyć tą średnią) dla anemometru CH4 powyższe równanie sprowadza się do zależności

\[ V_2=V(h_2) = V_1 \cdot \left( \frac{h_2}{h_1} \right)^w \quad (2) \]

i zawiera jedną niewiadomą \(w\), która jesteśmy wyznaczyć (jakiś logarytm o jakiejś podstawie się przyda, sprawdzić działanie polecenia log w pomocy).

  1. Wyprowadzić wzór na współczynnik \(w\) korzystając z równania (2) - “na papierze”

  2. Na podstawie danych pomiarowych (średnich) CH4 i CH6 wyznaczyć wartość współczynnika \(w\), (powinno wyjść z grubsza wartość z zakresu 0.2-0.3, tak tylko orientacyjnie).

  3. Zbudować funkcję jak we wzorze (1), argumentem powinna być tylko wysokość \(h\), a \(w\), \(h_1\) i \(V_1\) ustawić jako stałe.

  4. Wyznaczyć za pomocą funkcji teoretyczną średnią prędkość wiatru na wysokości odpowiadającej CH5 i porównać z danymi pomiarowymi CH5

  5. Jakiej, teoretycznie, wartości można się spodziewać na wysokości 130m (jest to często spotykana wysokość montażu gondoli turbiny wiatrowej - tam gdzie zaczepione są śmigła wiatraka).

  6. Narysować wykres tej funkcji dla \(h \in [10,130]\)

  7. Nanieść na wykres wartości z pomiaru (dla CH4 CH5 i CH6) - wskazówka: może się przydać konstrukcja wykresu z zadania 6.1 z Laboratorium 6

Wykres (wraz z częścią z następnego zadania) powinien przypominać

2.3 Zadanie - prawo logarytmiczne profilu wiatru

Inny sposób opisania zależności prędkości wiatru od wysokości to prawo logarytmiczne:

\[ V(h) = V_1 \cdot \frac{\ln(h/z_0)}{\ln(h_1/z_0)} \]

gdzie:

\(h_1\) - bazowa wysokość
\(V_1\) - średnia prędkość zmierzona na bazowej wysokości
\(z_0\) - współczynnik szorstkości terenu (roughness length). Jest on charakterystyczny dla danego miejsca pomiaru, podobnie jak wcześniej można go wyliczyć mając wartości pomiaru dla dwóch wysokości. Ma wartości zazwyczaj w przedziale 0.0002 dla obszarów wodnych (jeziora, morze) do 1.6 przy gęstej zabudowie.

Zadanie: ponowić obliczenia jak w poprzednim zadaniu, tym razem przy wykorzystaniu prawa logarytmicznego. \(z_0 \approx 1.3\) - tak mniej więcej powinno wyjść. Wykres prawa logarytmicznego nanieść na wykres z prawa Hellmanna.

Więcej informacji:

https://www.intechopen.com/chapters/17121

2.4 Wzorki, jakby ktoś nie umiał wyprowadzić:

Mając znane dwie wysokości \(h_1, h_2\) i policzone średnie prędkości \(V_1, V_2\) możemy wyliczyć

\[ w = \log_{\frac{h_2}{h_1}}{\frac{V_2}{V_1}} \]

albo równoważnie korzystając z własności logarytmu:

\[ w= \frac{\ln V_2-\ln V_1}{\ln h_2-\ln h_1} \]

oraz

\[ z_0=\left(\frac{h_1^{V_2/V_1}}{h_2}\right)^{\frac{1}{V_2/V_1-1}} \]

albo

\[ z_0=\exp\left(\frac{h_1^w \cdot \ln h_2-h_2^w \cdot \ln h_1}{h_1^w-h_2^w}\right) \]

3 Profil poziomy wiatru

3.1 Zadanie - różnica wskazań vane’ów

Instalacje nie są idealne, każde “skrzywienie” może wpłynąć na pomiar kierunku wiatru.

  • Narysować wykres wskazań vane’ów (wiatrowskazów) CH7 i CH8 w zależności od czasu. Powinna być widoczna różnica we wskazaniach.

  • Narysować wykres wartości bezwzględnej różnicy wskazań vane’ów CH7 i CH8. Wskazówka: można obliczyć tą różnicę i zapisać jako dodatkowy element ramki jak w Laboratorium 7, pakiet dyplr, zadanie z powierzchnią płatka

  • Ocenić średnią różnicę między wskazaniami vane’ów, biorąc pod uwagę że mamy do czynienia z kątami z zakresu od 0 do 360 stopni. Przykładowo, kierunki mające wartość 359 i 1st. są właściwie takie sam, choć ich różnica jest bardzo duża. Dodatkowe utrudnienie: przy niskich prędkościach wiatru kierunek może być bardzo zmienny w czasie i między wysokościami (kto puszczał latawce, ten wie), więc wypadałoby wykluczyć z analizy prędkości poniżej 3 m/s dla CH6. Ponownie tutaj można wykorzystać pakiet dyplr, gdzie funkcja mutate w połączeniu case_when w miejsce niezdefiniowanego wyniku wstawia NA

  • Narysować ponownie wykres wskazań vane’ów (wiatrowskazów) CH7 i CH8 w zależności od czasu, tym razem z uwzględnieniem wyliczonej różnicy (np. rozdzielając ją po połowie i dodając do CH8Avg i odejmując od CH7Avg tą część różnicy).

3.1.1 Przykładowy wynik:

## [1] "Obliczona poprawka to (przykladowo)"
## [1] 19.12772

3.2 Dopasowanie rozkładu do danych

Istotnym elementem analizy statystycznej danych jest dopasowanie rozkładu do danych empirycznych. Ograniczymy się tutaj do gotowej procedury z pakietu fitdistrplus.

library(fitdistrplus)

Dla przykładu wygenerujemy 1000 obserwacji z rozkładu normalnego i sprawdzimy jak działa funkcja dopasowująca rozkład dla trzech rodzin rozkładów (poprawny rozkład normalny, dodatkowe rozkłady jednostajny i gamma dla porównania):

testowy = rnorm(1000, mean=8, sd=2)
dopasowanie1 = fitdist(testowy, "norm") #podajemy wektor testowy oraz rodzinę rozkładów, jaka nas interesuje
dopasowanie2 = fitdist(testowy, "unif") #wynikiem jest obiekt klasy fitdist
dopasowanie3 = fitdist(testowy, "gamma")
dopasowanie1
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters:
##      estimate Std. Error
## mean 7.957226 0.06232158
## sd   1.970781 0.04406796
dopasowanie2
## Fitting of the distribution ' unif ' by maximum likelihood 
## Parameters:
##     estimate
## min  1.02928
## max 14.32074
dopasowanie3
## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters:
##        estimate Std. Error
## shape 14.390615  0.6362445
## rate   1.808566  0.0813700

Funkcja gofstat pozwala na wyświetlenie oceny dopasowania parametrów poprzez różne kryteria. Generalnie, im mniejsza wartość tym lepsze dopasowanie:

gofstat(list(dopasowanie1, dopasowanie2, dopasowanie3)) #ocena dopasowania poprzez różne kryteria - ogólnie, im mniej tym lepsze dopasowanie
## Goodness-of-fit statistics
##                              1-mle-norm 2-mle-unif 3-mle-gamma
## Kolmogorov-Smirnov statistic 0.03410812  0.2308232  0.06658204
## Cramer-von Mises statistic   0.20807755 22.0364427  1.16008348
## Anderson-Darling statistic   1.09835439        Inf  6.11363741
## 
## Goodness-of-fit criteria
##                                1-mle-norm 2-mle-unif 3-mle-gamma
## Akaike's Information Criterion   4198.737   5178.244    4276.221
## Bayesian Information Criterion   4208.553   5188.060    4286.037

Sposób na wyłuskanie z wyników samych parametrów:

dopasowanie1$estimate[1] # sposób na wyświetlenie parametru
##     mean 
## 7.957226
dopasowanie1$estimate[2] # sposób na wyświetlenie parametru
##       sd 
## 1.970781

W ten sposób wyłuskane parametry można użyć bezpośrednio lub przypisać do zmiennej pomocniczej.

3.3 Zadanie - histogram i rozkład prawdopodobieństwa

Do modelowania rozkładu prędkości wiatru najczęściej stosuje się rozkład Weibulla - informacje o rozkładzie Weibulla, który ma dwa parametry: kształt shape i skala scale.

Zajmiemy się analizą prędkości wiatru dla anemometru CH4.

  1. Sprawdzić podstawowe informacje dotyczące CH4 poprzez polecenie summary

  2. Narysować histogram prędkości wiatru z podziałem co 0.5 m/s. Może się przydać parametr binwidth. Będziemy dopasowywać gęstość i będziemy chcieli nanieść ja na histogram, zatem może przydać się zamiana liczności na częstotliwość występowania poprzez użycie y = after_stat(count / sum(count)) w funkcji aes histogramu.

  3. Za pomocą pakietu fitdistrplus wyestymować parametry rozkładu weibull.

  4. Nanieść wykres gęstości z otrzymanymi parametrami na wykres histogramu. Można wykorzystać geometrię: geom_function(fun = function(x) (dweibull(x, shape = ksztalt, scale = skala)), linewidth = 1, colour = "black") gdzie w miejsce ksztalt i skala należy podać wyestymowane wartości

Wykres powinien przypominać:

UWAGA: możliwe jest też przeskalowanie samej funkcji rozkładu Weibulla poprzez przemnożenie dweibull(x, shape = ksztalt, scale = skala) przez liczbę obserwacji razy szerokość okna histogramu:

UWAGA: druga skala z procentami została uzyskana przez scale_y_continuous(sec.axis = sec_axis(~./(4464), name = "Gęstość rozkładu",labels = scales::percent), name = "Liczba wystąpień")

3.4 Zadanie - parametry rozkładu Weibulla - cd

  1. Załadować bibliotekę library("mixdist") . Może być potrzebne jej jednorazowe zainstalowanie.

  2. Powyższą bibliotekę wykorzystamy do wyznaczenia średniej (i odchylenia standardowego)
    weibullparinv(shape=???, scale=???, loc = 0),
    więcej na
    https://search.r-project.org/CRAN/refmans/mixdist/html/weibullparinv.html

  3. Czy uzyskana jak wyżej średnia zgadza się ze średnią empiryczną? (Średnia empiryczna - zwykła średnia dla CH4)

  4. Cel tej analizy to oszacowanie, ile energii można wyprodukować w danym miejscu.
    Pod linkiem

    Przystępny kalkulator

    dostępny jest kalkulator.

    W oknie kalkulatora należy ustawić parametr \(A\) - nasz parametr skali i \(k\) - nasz parametr kształtu.


    Poeksperymentować z różnymi modelami turbin.

3.4.1 Gęstość powietrza

Gęstość powietrza można pozostawić bez zmian wyliczyć na podstawie danych ciśnienia - CH10 w hPa oraz temperatury CH11 w stopniach Celcjusza, korzystając ze wzoru: \[ \rho=\frac{p \cdot M}{R \cdot T} \]

gdzie:

p - ciśnienie w Pa (trzeba przeliczyć, 1hPa=100Pa)
M - efektywna masa molowa, 0.0289 kg/mol dla powietrza
R - uniwersalna stała gazowa równa 8.31446261815324 J/(mol*K)
T - temperatura w skali bezwzględnej (Kelvinach), trzeba przeliczyć
Powinno wyjść około 1.2 kg/m3