MATLAB Examples

Description

In this script, I reproduce the results presented by John D. Holmes in the first part of the chapter 2 of his book: Wind loading of structures [1]. The notations he uses are slightly different in his book than in the Matlab statistics toolbox. In particular, the expression of the cumulative density function of the GEV he uses is based on the shape factor k whereas the statistical toolbox uses "K = -k". To keep consistent with the Matlab notations, I will also use the convention "K = -k" in the following.

[1] Holmes, J. D. (2015). Wind loading of structures. CRC Press.

Contents

Introduction

The figure plotted here reproduces the one shown in Fig 2.1 in [1]

clearvars;close all;clc;
K = [0,0.2,-0.2]; % 3 types of GEV distributions (K<0, K=0 and K>0)
a = 1;
u = 1;
U = linspace(0,2,100);

% Check if statistic toolbox exists
status = license('test','statistics_toolbox');
if status==0, warning('No satistics toolbox detected');end

clear F*
for ii=1:numel(K),
    F(ii,:) = mygevcdf(U,K(ii),a,u); % if no stat toolbox available
end

subplot(211)
hold on;box on;
plot(-log(-log(F(1,:))),(U-u)./a,'k','linewidth',2);
plot(-log(-log(F(2,:))),(U-u)./a,'r','linewidth',0.7);
plot(-log(-log(F(3,:))),(U-u)./a,'b:','linewidth',2);
grid on
xlabel('-ln(-ln(Fu))')
    ylabel('(U-u)/a')
legend('Type I k = 0','Type II k = -0.2','Type III k = +0.2',...
    'location','Northwest')
title('no stat toolbox available: use of mygevcdf')

if status ==1,
    for ii=1:3,
        F(ii,:)  = gevcdf(U,K(ii),a,u); % if stat toolbox is available
    end

    subplot(212)
    hold on;box on;
    plot(-log(-log(F(1,:))),(U-u)./a,'k','linewidth',2);
    plot(-log(-log(F(2,:))),(U-u)./a,'r','linewidth',0.7);
    plot(-log(-log(F(3,:))),(U-u)./a,'b:','linewidth',2);
    grid on
    xlabel('-ln(-ln(Fu))')
    ylabel('(U-u)/a')
    title(' stat toolbox available: use of gevcdf')
end

set(gcf,'color','w')

Example 1

The data used (EastSale.mat) in this example comes from page 42-43 in [1]. The extreme 10 min win velocities are calculated for different return periods, and different methods. The matlab built-in function provided with the statistical toolbox is compared to the one used here.

load EastSale
if status ==1,
    [parmhat0,parmci0] = gevfit(meanU); % Matlab built-in function
else
    parmhat0 = nan;
    parmci0 = nan;
end

% Home-made function
[parmhat1] = fitGEV(meanU,'method','Gumbel');
[parmhat2] = fitGEV(meanU,'method','moments');
[parmhat3] = fitGEV(meanU,'method','Gringorten');

fprintf('Parameters found: \n')
disp(['--------------------------------------------------------------------------------------------'])
disp(['Method          ','Gumbel' ,'     moments','     Gringorten','     gevfit (Matlab)'])
disp(['--------------------------------------------------------------------------------------------'])
disp(['k','               ',num2str(parmhat1(1),4),'         ',num2str(parmhat2(1),4),'            ',num2str(parmhat3(1),4),'             ',num2str(parmhat0(1),4)])
disp(['mu','              ',num2str(parmhat1(2),4),'     ',num2str(parmhat2(2),4),'        ',num2str(parmhat3(2),4),'          ',num2str(parmhat0(2),4)])
disp(['sigma','           ',num2str(parmhat1(3),4),'     ',num2str(parmhat2(3),4),'        ',num2str(parmhat3(3),4),'          ',num2str(parmhat0(3),4)])
disp(['--------------------------------------------------------------------------------------------'])
fprintf('\n\n')
% cdf of GEV
X = linspace(min(meanU),3*max(meanU),1e4); % Measured
F0  = mygevcdf(X,parmhat0(1),parmhat0(2),parmhat0(3)); % matlab
F1  = mygevcdf(X,parmhat1(1),parmhat1(2),parmhat1(3)); % Gumbel
F2  = mygevcdf(X,parmhat2(1),parmhat2(2),parmhat2(3)); % moments
F3  = mygevcdf(X,parmhat3(1),parmhat3(2),parmhat3(3)); % Gringorten
% Return period R
R(1,:)= 1./(1-F1); % Gumbel
R(2,:) = 1./(1-F3); % Gringorten
R(3,:) = 1./(1-F2); % moments
R(4,:) = 1./(1-F0); % matlab
clear ind*
% Find return period
for ii=1:4
    [~,ind_50(ii)] = min(abs(R(ii,:)-50));
    [~,ind_100(ii)] = min(abs(R(ii,:)-100));
    [~,ind_200(ii)] = min(abs(R(ii,:)-200));
    [~,ind_500(ii)] = min(abs(R(ii,:)-500));
    [~,ind_1000(ii)] = min(abs(R(ii,:)-1000));
end

fprintf('Prediction of extreme wind speeds for East Sale (synoptic winds): \n')


val = [X(ind_50);X(ind_100);X(ind_200);X(ind_500);X(ind_1000)];
disp(['--------------------------------------------------------------------------------------------'])
disp(['Return period (years)','              Predicted       gust       speed       (m/s)'])
disp(['--------------------- ',' ------------------------------------------------------------------'])
disp(['       ','                  (Gumbel)','       (Gringorten)','       (moments)','       (Matlab)'])
disp(['50','                        ',num2str(val(1,1),3),'           ',num2str(val(1,2),3),'               ',num2str(val(1,3),3),'             ',num2str(val(1,4),3)])
disp(['100','                       ',num2str(val(2,1),3),'             ',num2str(val(2,2),3),'               ',num2str(val(2,3),3),'             ',num2str(val(2,4),3)])
disp(['200','                       ',num2str(val(3,1),3),'           ',num2str(val(3,2),3),'                 ',num2str(val(3,3),3),'             ',num2str(val(3,4),3)])
disp(['500','                       ',num2str(val(4,1),3),'           ',num2str(val(4,2),3),'               ',num2str(val(4,3),3),'             ',num2str(val(4,4),3)])
disp(['1000','                      ',num2str(val(5,1),3),'           ',num2str(val(5,2),3),'                 ',num2str(val(5,3),3),'             ',num2str(val(5,4),3)])
disp(['--------------------------------------------------------------------------------------------'])
Parameters found: 
--------------------------------------------------------------------------------------------
Method          Gumbel     moments     Gringorten     gevfit (Matlab)
--------------------------------------------------------------------------------------------
k               0         0            0             -0.001661
mu              2.659     2.492        2.513          2.421
sigma           27.81     27.83        27.84          27.89
--------------------------------------------------------------------------------------------


Prediction of extreme wind speeds for East Sale (synoptic winds): 
--------------------------------------------------------------------------------------------
Return period (years)              Predicted       gust       speed       (m/s)
---------------------  ------------------------------------------------------------------
                         (Gumbel)       (Gringorten)       (moments)       (Matlab)
50                        38.2           37.6               37.6             37.3
100                       40             39.4               39.3             39
200                       41.9           41.1                 41             40.7
500                       44.3           43.4               43.3             42.9
1000                      46.2           45.2                 45             44.5
--------------------------------------------------------------------------------------------

Data plot

This time the function mygevcdf is used in combination with fitGEV.
% get cdf
returnPeriod = 50;
x = ([1:numel(meanU)])/(numel(meanU)+1);
r = 1./(1-x);
X = linspace(min(meanU),2*max(meanU),1e4);

F  = mygevcdf(X,parmhat2(1),parmhat2(2),parmhat2(3)); % for fitted parameters
R = 1./(1-F);

figure
hold on;box on;
plot(r,sort(meanU),'ko','markerfacecolor','r'); % cdf of meanU
plot(R,X,'k');
xlabel('return period (years)');
ylabel('wind velocity (m/s)');
axis tight
  ylim([0.9*min(X(R<returnPeriod)),1.1*max(X(R<returnPeriod))])
    xlim([0,returnPeriod])
legend('Measured','fitted','location','SouthEast')
grid on
set(gcf,'color','w');

More straightforward data plot

It is actually possible to plot directly the figure using fitGEV:
load EastSale
[parmhat2] = fitGEV(meanU,'method','moments','dataPlot',1,'returnPeriod',50);
50-year hourly wind velocity is 37.6 m/s (based on 56 years of data)