% This function loads the variable 'MaxEvents' which is produced by 'MATSimAxles'
load('Q1Q2MaxEventsx.mat');
% Convert ClassT to m (number) form m = 1 is All, m = 2 is ClassOW, m = 3 is Class
MaxEvents.m(strcmp(MaxEvents.ClassT,"All")) = 1;
MaxEvents.m(strcmp(MaxEvents.ClassT,"ClassOW")) = 2;
MaxEvents.m(strcmp(MaxEvents.ClassT,"Class")) = 3; MaxEvents.ClassT = [];
MaxEvents.L1Sp = MaxEvents.L1Sp/100;
MaxEvents.L2Sp = MaxEvents.L2Sp/100;
% Transform AxTandem into Array (necessary for splitapply)
Z.Time = datenum(Z.Time);
for i = 1:length(ClassType)
% Filter based on Class - MaxEvents is compromised after this (deleting entries)
if strcmp(Class,'ClassOW')
MaxEvents(MaxEvents.m == 1,:) = []; %#ok<*SAGROW>
elseif strcmp(Class,'Class')
MaxEvents(MaxEvents.m == 1,:) = [];
MaxEvents(MaxEvents.m == 2,:) = [];
Max.(Class).(BlockM) = [];
if strcmp(BlockM,'Daily')
% Make groups out of unique locations and days
[Gr, GrIDDay, GrIDZST] = findgroups(dateshift(MaxEvents.Time,'start','day'),MaxEvents.ZST);
elseif strcmp(BlockM,'Weekly')
[Gr, GrIDWeek, GrIDZST] = findgroups(dateshift(MaxEvents.Time,'start','week'),MaxEvents.ZST);
[Gr, GrIDYear, GrIDZST] = findgroups(year(MaxEvents.Time),MaxEvents.ZST);
% Perform splitapply (see function at end... not just Max as we want whole rows involving maxes)
Max.(Class).(BlockM) = splitapply(@(Z)maxIndex(Z),Z,Gr);
% Transform back into table form
Max.(Class).(BlockM) = array2table(Max.(Class).(BlockM));
Max.(Class).(BlockM).Properties.VariableNames = {'R', 'ZST', 'Max', 'DayRank', 'L1Veh', 'L2Veh', 'L1Load', 'L2Load', 'L1Ax', 'L2Ax', 'L1Sp', 'L2Sp', 'Time', 'm'};
Max.(Class).(BlockM).Time = datetime(Max.(Class).(BlockM).Time,'ConvertFrom',"datenum"); Max.(Class).(BlockM).R = [];
DistTypes = {'Normal', 'Lognormal', 'GeneralizedExtremeValue'};
% Set CDF Scaling Factors for estimates
D2WFactor = 5; W2YFactor = 50; D2YFactor = D2WFactor*W2YFactor;
% Set x values to be use globally
% Fit Block Maxima to Normal Curve
for i = 1:length(ClassType)
for k = 1:length(DistTypes)
pd.(Class).(BlockM).(Dist) = fitdist(Max.(Class).(BlockM).Max,Dist);
y_values.(Class).(BlockM).(Dist).PDF_Fit = pdf(pd.(Class).(BlockM).(Dist),x_values);
y_values.(Class).(BlockM).(Dist).CDF_Fit = cdf(pd.(Class).(BlockM).(Dist),x_values);
if strcmp(BlockM,'Daily')
y_values.(Class).('Weekly').(Dist).CDF_D2W = y_values.(Class).(BlockM).(Dist).CDF_Fit.^D2WFactor;
y_values.(Class).('Weekly').(Dist).PDF_D2W = [0 diff(y_values.(Class).('Weekly').(Dist).CDF_D2W)];
y_values.(Class).('Yearly').(Dist).CDF_D2Y = y_values.(Class).(BlockM).(Dist).CDF_Fit.^D2YFactor;
y_values.(Class).('Yearly').(Dist).PDF_D2Y = [0 diff(y_values.(Class).('Yearly').(Dist).CDF_D2Y)];
mu = trapz(y_values.(Class).('Weekly').(Dist).CDF_D2W,x_values);
sigma = sqrt(trapz(y_values.(Class).('Weekly').(Dist).CDF_D2W,(x_values-mu).^2));
pd.(Class).D2W.(Dist) = makedist(Dist,"mu",mu,"sigma",sigma);
mu = trapz(y_values.(Class).('Yearly').(Dist).CDF_D2Y,x_values);
sigma = sqrt(trapz(y_values.(Class).('Yearly').(Dist).CDF_D2Y,(x_values-mu).^2));
pd.(Class).D2Y.(Dist) = makedist(Dist,"mu",mu,"sigma",sigma);
elseif strcmp(Dist,'Lognormal')
mu = trapz(y_values.(Class).('Weekly').(Dist).CDF_D2W,x_values);
sigma = sqrt(trapz(y_values.(Class).('Weekly').(Dist).CDF_D2W,(x_values-mu).^2));
zeta = sqrt(log(1+(sigma/mu)^2));
lambda = log(mu)-0.5*zeta^2;
pd.(Class).D2W.(Dist) = makedist(Dist,"mu",lambda,"sigma",zeta);
mu = trapz(y_values.(Class).('Yearly').(Dist).CDF_D2Y,x_values);
sigma = sqrt(trapz(y_values.(Class).('Yearly').(Dist).CDF_D2Y,(x_values-mu).^2));
zeta = sqrt(log(1+(sigma/mu)^2));
lambda = log(mu)-0.5*zeta^2;
pd.(Class).D2Y.(Dist) = makedist(Dist,"mu",lambda,"sigma",zeta);
elseif strcmp(BlockM,'Weekly')
y_values.(Class).('Yearly').(Dist).CDF_W2Y = y_values.(Class).(BlockM).(Dist).CDF_Fit.^W2YFactor;
y_values.(Class).('Yearly').(Dist).PDF_W2Y = [0 diff(y_values.(Class).('Yearly').(Dist).CDF_W2Y)];
mu = trapz(y_values.(Class).('Yearly').(Dist).CDF_W2Y,x_values);
sigma = sqrt(trapz(y_values.(Class).('Yearly').(Dist).CDF_W2Y,(x_values-mu).^2));
pd.(Class).W2Y.(Dist) = makedist(Dist,"mu",mu,"sigma",sigma);
elseif strcmp(Dist,'Lognormal')
mu = trapz(y_values.(Class).('Yearly').(Dist).CDF_W2Y,x_values);
sigma = sqrt(trapz(y_values.(Class).('Yearly').(Dist).CDF_W2Y,(x_values-mu).^2));
zeta = sqrt(log(1+(sigma/mu)^2));
lambda = log(mu)-0.5*zeta^2;
pd.(Class).W2Y.(Dist) = makedist(Dist,"mu",lambda,"sigma",zeta);
figure('Position',[0 0 1500 400]);
x = X(1:end-1) + diff(X);
for i = 1:length(ClassType)
y = histcounts(Max.(Class).(BlockM).Max,'BinEdges',X,'normalization','pdf');
bar(x,y/ScaleDown(j),1,'EdgeColor',C(j,:),'FaceColor',[.8 .8 .8],'FaceAlpha',0.5,'DisplayName',BlockM)
plot(x_values,y_values.(Class).(BlockM).(Dist).PDF_Fit/ScaleDown(j),'k-','DisplayName',[Dist 'Fit'])
if strcmp(BlockM,'Yearly')
plot(x_values,y_values.(Class).(BlockM).(Dist).PDF_D2Y/ScaleDown(j),'k--','DisplayName',[Dist 'D2Y'])
plot(x_values,y_values.(Class).(BlockM).(Dist).PDF_W2Y/ScaleDown(j),'k:','DisplayName',[Dist 'W2Y'])
set(gca,'ytick',[],'yticklabel',[],'ycolor','k')
xlim([LimitL+25 LimitR-25])
ylabel('Normalized Histograms')
xlabel('Total Strip Load (kN)')
title(['Maximum Strip Load ' Class])
legend('location','northeast')
figure('Position',[0 0 1500 400]);
Xp = LimitL-Step/2:Step:LimitR-Step/2;
histogram(Max.All.(BlockM).Max,'BinEdges',Xp,'normalization','pdf','EdgeColor',C(4,:),'FaceColor',[.8 .8 .8],'FaceAlpha',0.2,'DisplayName','All')
histogram(Max.ClassOW.(BlockM).Max,'BinEdges',Xp,'normalization','pdf','EdgeColor',C(7,:),'FaceColor',[.8 .8 .8],'FaceAlpha',0.2,'DisplayName','ClassOW')
histogram(Max.Class.(BlockM).Max,'BinEdges',Xp,'normalization','pdf','EdgeColor',C(9,:),'FaceColor',[.8 .8 .8],'FaceAlpha',0.2,'DisplayName','Class')
set(gca,'ytick',[],'yticklabel',[])
ylabel('Normalized Histograms')
xlabel('Total Strip Load (kN)')
title(['Maximum Tandem Axles ' BlockM])
legend('location','best')
xlim([LimitL+25 LimitR-25])
%sum(Max.All.Daily.L1Load)./sum(Max.All.Daily.L2Load)
%histogram(hour(Max.All.Daily.Time),'normalization','pdf')
% All TrTyps in number format
TrTyps = [11; 12; 22; 23; 111; 11117; 1127; 12117; 122; 11127; 1128; 1138; 1238; 41; 42; 43; 44; 45; 46; 0; 99];
TrName{1} = '11'; TrName{2} = '12'; TrName{3} = '22'; TrName{4} = '23'; TrName{5} = '111'; TrName{6} = '1111r'; TrName{7} = '112r';
TrName{8} = '1211r'; TrName{9} = '122'; TrName{10} = '1112r'; TrName{11} = '112a'; TrName{12} = '113a'; TrName{13} = '123a';
TrName{14} = '5ax 60t'; TrName{15} = '6ax 60t'; TrName{16} = '7ax 72t'; TrName{17} = '8ax 84t'; TrName{18} = '9ax 96t';
TrName{19} = '8ax 96t'; TrName{20} = 'No Class'; TrName{21} = 'Empty'; TrName = TrName';
% Convert from Cell to String
TN = strings(size(TrName)); [TN{:}] = TrName{:};
% For each ClassType (or in this case just All and Class)
% Select Daily, Weekly, or Yearly
figure('Position',[0 0 1500 400]);
% Get Primary Lane and Accompanying Lane Vehicles
Primary = zeros(height(Max.(Class).(BlockM)),1); Accomp = Primary;
[~, Ind] = max([Max.(Class).(BlockM).L1Load'; Max.(Class).(BlockM).L2Load']); Ind = Ind';
Primary(Ind == 1) = Max.(Class).(BlockM).L1Veh(Ind == 1); Primary(Ind == 2) = Max.(Class).(BlockM).L2Veh(Ind == 2);
Accomp(Ind == 1) = Max.(Class).(BlockM).L2Veh(Ind == 1); Accomp(Ind == 2) = Max.(Class).(BlockM).L1Veh(Ind == 2);
PrimaryC = categorical(Primary); AccompC = categorical(Accomp);
[N, Cats] = histcounts(PrimaryC,categorical(TrTyps)); [Nx] = histcounts(AccompC,categorical(TrTyps));
Label = cell(length(TN),1); Labelx = cell(length(TN),1);
Perc = 100*N/sum(N); Percx = 100*Nx/sum(Nx);
PercS = string(round(Perc,1))'; PercSx = string(round(Percx,1))';
Label{p} = [TN{p} ' - ' PercS{p} '%']; else; Label{p} = '';
Labelx{p} = [TN{p} ' - ' PercSx{p} '%']; else; Labelx{p} = '';
% We need colors to be the same for each TrTyp
% Labels should be strings that include the percentages
sgtitle([Class ' ' BlockM],'fontweight','bold','fontsize',12);
newColors = linspecer(numel(N));
% Isolate the patch handles
patchHand = findobj(h, 'Type', 'Patch');
% Set the color of all patches using the nx3 newColors matrix
set(patchHand, {'FaceColor'}, mat2cell(newColors, ones(size(newColors,1),1), 3))
% Isolate the patch handles
patchHand = findobj(h, 'Type', 'Patch');
set(patchHand, {'FaceColor'}, mat2cell(newColors, ones(size(newColors,1),1), 3))
title('Accompanying Lane')
if strcmp(Class,'All') && strcmp(BlockM,'Daily')
NAD = N + Nx; NAD = NAD/sum(NAD);
elseif strcmp(Class,'All') && strcmp(BlockM,'Weekly')
NAY = N + Nx; NAY = NAY/sum(NAY);
% One question is: how much more prevalent is the TrTyp in the maximums
% than in the normal traffic? Who are the heavy ppl? Use ClassOW
% Need to fetch avg prevalance of TrTyps from somewhere...
% Lets do the comparison between AllDaily and AllYearly
figure('Position',[0 0 1500 400]);
b = bar(categorical(TN),100*NAY./NAD);
b.CData(:,:) = newColors;
ylabel('Prevalance Increase (%)')
title('Prevalance Increase (%) From Daily to Weekly')
b = bar(categorical(TN),100*NAY./NAD);
b.CData(:,:) = newColors;
ylabel('Prevalance Increase (%)')
title('Prevalance Increase (%) From Daily to Weekly')
figure('Position',[0 0 1500 400]);
%Xp = LimitL-Step/2:Step:LimitR-Step/2;
x = X(1:end-1) + diff(X);
for i = 1:length(ClassType)
y = histcounts(Max.(Class).(BlockM).L1Sp(Max.(Class).(BlockM).L1Sp>0),'BinEdges',X,'normalization','pdf');
bar(x,y/ScaleDown(j),1,'EdgeColor',C(j,:),'FaceColor',[.8 .8 .8],'FaceAlpha',0.5,'DisplayName',BlockM)
set(gca,'ytick',[],'yticklabel',[],'ycolor','k')
xlim([LimitL+.25 LimitR-.25])
ylabel('Normalized Histograms')
title(['Maximum Strip Vehs - Speed ' Class])
legend('location','northwest')
% According to Annex C of SIA 269
Beta.Yearly = 4.7; % Here we use the Beta annual - so we should use annual max effects
Beta.Weekly = 5.44422; % See Tail Fitting > Beta Conversion
% Right now we assume a normal distribution in the Edact calculation. Also
SumTableNorm = array2table(zeros(6,length(BM)*length(ClassType)));
SumTableLogNormComp = array2table(zeros(5,length(BM)*length(ClassType)));
SumTableNorm.Properties.VariableNames = {[BM{3} ' ' ClassType{3}], [BM{3} ' ' ClassType{2}], [BM{3} ' ' ClassType{1}], [BM{2} ' ' ClassType{3}],...
[BM{2} ' ' ClassType{2}], [BM{2} ' ' ClassType{1}], [BM{1} ' ' ClassType{3}], [BM{1} ' ' ClassType{2}], [BM{1} ' ' ClassType{1}]};
SumTableLogNorm = SumTableNorm;
SumTableNorm.Properties.RowNames = {'Mu', 'Sigma', 'COV (%)', '95th %', 'Dist 95', 'AlphaQ'};
SumTableLogNorm.Properties.RowNames = {'Lambda', 'Zeta', '95th %', 'Dist 95', 'AlphaQN', 'AlphaQLN'};
SumTableLogNormComp.Properties.RowNames = {'Lambda', 'Zeta', '95th %', 'Dist 95', 'AlphaQLN Est'};
SumTableLogNormComp.Properties.VariableNames = {[BM{3} ' ' ClassType{3}], [BM{3} ' ' ClassType{2}], [BM{3} ' ' ClassType{1}], ['W2Y ' ClassType{3}],...
['W2Y ' ClassType{2}], ['W2Y ' ClassType{1}], ['D2Y ' ClassType{3}], ['D2Y ' ClassType{2}], ['D2Y ' ClassType{1}]};
for i = 1:length(ClassType)
Emact = mean(Max.(Class).(BlockM).Max);
Stdev = std(Max.(Class).(BlockM).Max);
N5 = prctile(Max.(Class).(BlockM).Max,95);
Dist95N = norminv(0.95,pd.(Class).(BlockM).Normal.mu,pd.(Class).(BlockM).Normal.sigma);
Dist95LN = logninv(0.95,pd.(Class).(BlockM).Lognormal.mu,pd.(Class).(BlockM).Lognormal.sigma);
EdactN = Emact*(1+Alpha*Beta.(BlockM)*COV);
AlphaQN = EdactN/(1000*1.5);
EdactLN = Emact*exp(Alpha*Beta.(BlockM)*sqrt(Delta2)-0.5*Delta2);
AlphaQLN = EdactLN/(1000*1.5);
SumTableNorm.([BlockM ' ' Class]) = [Emact; Stdev; 100*COV; N5; Dist95N; AlphaQN];
SumTableLogNorm.([BlockM ' ' Class]) = [pd.(Class).(BlockM).Lognormal.mu; pd.(Class).(BlockM).Lognormal.sigma; N5; Dist95LN; AlphaQN; AlphaQLN];
if strcmp(BlockM,'Weekly')
pdx = pd.(Class).W2Y.Lognormal;
Emact = exp(pdx.mu+0.5*Delta2);
SumTableLogNormComp.(['W2Y ' Class]) = [pdx.mu; pdx.sigma; prctile(Max.(Class).Yearly.Max,95); logninv(0.95,pdx.mu,pdx.sigma); Emact*exp(Alpha*Beta.Yearly*sqrt(Delta2)-0.5*Delta2)/(1000*1.5)];
elseif strcmp(BlockM,'Daily')
pdx = pd.(Class).D2Y.Lognormal;
Emact = exp(pdx.mu+0.5*Delta2);
SumTableLogNormComp.(['D2Y ' Class]) = [pdx.mu; pdx.sigma; prctile(Max.(Class).Yearly.Max,95); logninv(0.95,pdx.mu,pdx.sigma); Emact*exp(Alpha*Beta.Yearly*sqrt(Delta2)-0.5*Delta2)/(1000*1.5)];
SumTableLogNormComp.([BlockM ' ' Class]) = [pd.(Class).(BlockM).Lognormal.mu; pd.(Class).(BlockM).Lognormal.sigma; N5; Dist95LN; AlphaQLN];
disp(SumTableNorm)
Yearly Class Yearly ClassOW Yearly All Weekly Class Weekly ClassOW Weekly All Daily Class Daily ClassOW Daily All
____________ ______________ __________ ____________ ______________ __________ ___________ _____________ _________
Mu 369.70 373.83 408.50 274.18 275.59 291.27 231.17 231.81 239.87
Sigma 48.93 50.51 68.04 36.69 37.72 44.97 35.54 36.02 40.79
COV (%) 13.23 13.51 16.66 13.38 13.69 15.44 15.37 15.54 17.00
95th % 444.04 460.94 520.67 344.25 347.88 378.33 289.40 290.38 308.43
Dist 95 450.18 456.91 520.41 334.53 337.64 365.24 289.62 291.06 306.96
AlphaQ 0.35 0.36 0.42 0.28 0.28 0.31 0.25 0.25 0.27
disp(SumTableLogNorm)
Yearly Class Yearly ClassOW Yearly All Weekly Class Weekly ClassOW Weekly All Daily Class Daily ClassOW Daily All
____________ ______________ __________ ____________ ______________ __________ ___________ _____________ _________
Lambda 5.90 5.91 6.00 5.61 5.61 5.66 5.43 5.43 5.47
Zeta 0.13 0.14 0.16 0.13 0.13 0.14 0.16 0.16 0.17
95th % 444.04 460.94 520.67 344.25 347.88 378.33 289.40 290.38 308.43
Dist 95 455.97 462.65 522.69 335.36 338.23 365.72 294.98 296.35 311.41
AlphaQN 0.35 0.36 0.42 0.28 0.28 0.31 0.25 0.25 0.27
AlphaQLN 0.38 0.38 0.46 0.30 0.31 0.34 0.28 0.28 0.31
disp(SumTableLogNormComp)
Yearly Class Yearly ClassOW Yearly All W2Y Class W2Y ClassOW W2Y All D2Y Class D2Y ClassOW D2Y All
____________ ______________ __________ _________ ___________ _______ _________ ___________ _______
Lambda 5.90 5.91 6.00 5.89 5.90 5.99 5.87 5.88 5.94
Zeta 0.13 0.14 0.16 0.06 0.06 0.07 0.06 0.06 0.07
95th % 444.04 460.94 520.67 444.04 460.94 520.67 444.04 460.94 520.67
Dist 95 455.97 462.65 522.69 399.99 404.70 446.98 392.31 395.07 423.28
AlphaQLN Est 0.38 0.38 0.46 0.29 0.30 0.33 0.29 0.29 0.32
figure('Position',[0 0 1500 400]);
for i = 1:length(ClassType)
[MaxECDF, MaxECDFRank] = ecdf(Max.(Class).(BlockM).Max); MaxECDFRank = MaxECDFRank'; MaxECDF(1) = []; MaxECDFRank(1) = [];
scatter(norminv((1:length(MaxECDFRank))/(length(MaxECDFRank) + 1)),log(MaxECDFRank),7,C(j,:),'filled','DisplayName','Max Data');
mdl = fitlm(norminv((1:length(MaxECDFRank))/(length(MaxECDFRank) + 1)),log(MaxECDFRank),'linear');
x2 = norminv((1:length(MaxECDFRank))/(length(MaxECDFRank) + 1));
md2 = fitlm(x2(round(length(x2)*3/4):length(x2)),y2(round(length(x2)*3/4):length(x2)),'linear');
plot(norminv((1:length(MaxECDFRank))/(length(MaxECDFRank) + 1)),mdl.Fitted,'-','Color',C(j,:),'DisplayName',['Fitted ' num2str(mdl.Rsquared.Ordinary,3)]);
%plot(x2(round(length(x2)*3/4):length(x2)),md2.Fitted,'-','Color','k','DisplayName',['Fitted ' num2str(md2.Rsquared.Ordinary,3)]);
%text(1,y1(1)+(y1(2)-y1(1))/j,['\lambda = ' sprintf('%.2f \n',mdl.Coefficients.Estimate(1)) '\zeta = ' sprintf('%.2f \n',mdl.Coefficients.Estimate(2)) sprintf('R^2: %.1f%%',mdl.Rsquared.Ordinary*100)])
text(1,y1(1)+(y1(2)-y1(1))*(j/6),sprintf('R^2: %.1f%%',mdl.Rsquared.Ordinary*100),"Color",C(j,:))
ylabel('log[ Total Strip Load (kN) ]')
xlabel('Standard Normal Percentile')
%legend('Location', 'southeast')
title(['Maximum Tandem Axles ' Class])
sgtitle('LogNormal Probability Paper','fontweight','bold','fontsize',12);