clear; clc; close all;
x=(0:1:150);
y=lognpdf(x,3.14,0.88);
m1=30.7;
v1=17.72;
numbmob=1681;
mu1 = log((m1^2)/sqrt(v1+m1^2));
sigma1 = sqrt(log(v1/(m1^2)+1));
d1=lognrnd(mu1,sigma1,numbmob,1);
m2=16;
v2=17.72;
numbmot=7743;
mu2 = log((m2^2)/sqrt(v2+m2^2));
sigma2 = sqrt(log(v2/(m2^2)+1));
d2=lognrnd(mu2,sigma2,numbmot,1);
Dmax1=150;
Eff1=0.95;
SOC1=(1-(d1/Dmax1*Eff1))*100;
Dmax2=50;
Eff2=0.95;
SOC2=(1-(d2/Dmax2*Eff2))*100;
CapMob=24;
ChargMob=6.6;
TC1=(0.95-(SOC1/100))*(CapMob/(ChargMob*Eff1));
CapMot=2;
ChargMot=0.55;
TC2=(0.95-(SOC2/100))*(CapMot/(ChargMot*Eff2));
for i=1:ceil(0.9*numbmob)
TPMob1(i,1)=normrnd(9,0.9);
while TPMob1(i,1)<9 || TPMob1(i,1)>17
TPMob1(i,1)=normrnd(9,0.9);
end
if TPMob1(i,1)>=24
TPMob1(i,1)=TPMob1(i,1)-24;
end
end
for i=1:numbmob-ceil(0.9*numbmob)
TPMob2(i,1)=normrnd(18.5,1);
while TPMob2(i,1)<18 || TPMob2(i,1)>31
TPMob2(i,1)=normrnd(18.5,1);
end
if TPMob2(i,1)>=24
TPMob2(i,1)=TPMob2(i,1)-24;
end
end
K1=[TPMob1;TPMob2];
TDMob=K1+TC1;
hh=0;
for h=0:0.5:23.5
hh=hh+1;
for i=1:length(K1)
if K1(i,1)<=h && TDMob(i,1)>=h
PEV1(i,hh)=2.2;
else
PEV1(i,hh)=0;
end
end
PEVT1(1,hh)=sum(PEV1(:,hh));
end
for i=1:ceil(0.9*numbmot)
TPMot1(i,1)=normrnd(9,0.9);
while TPMot1(i,1)<9 || TPMot1(i,1)>17
TPMot1(i,1)=normrnd(9,0.9);
end
if TPMot1(i,1)>=24
TPMot1(i,1)=TPMot1(i,1)-24;
end
end
for i=1:numbmot-ceil(0.9*numbmot)
TPMot2(i,1)=normrnd(18,5.1);
while TPMot2(i,1)<18 || TPMot2(i,1)>31
TPMot2(i,1)=normrnd(18,5.1);
end
if TPMot2(i,1)>=24
TPMot2(i,1)=TPMot2(i,1)-24;
end
end
K2=[TPMot1;TPMot2];
TDMot=K2+TC2;
hh=0;
for h=0:0.5:23.5
hh=hh+1;
for i=1:length(K2)
if K2(i,1)<=h && TDMot(i,1)>=h
PEV2(i,hh)=0.5;
else
PEV2(i,hh)=0;
end
end
PEVT2(1,hh)=sum(PEV2(:,hh));
end
PEVTotal=PEVT1+PEVT2;
T=([0:0.5:23.5]);
figure(1)
plot(x,y)
xticks([0:20:150])
title('Daily Travel Distance')
xlabel('Distance (km)')
ylabel('Probabilitas')
figure(2)
scatter(K1,SOC1,'*')
xticks([0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23])
yticks('auto')
title('Car Plug in Time')
xlabel('Hour')
ylabel('SOC')
figure(3)
scatter(K2,SOC2,'*')
xticks([0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23])
yticks('auto')
title('Motor Cycle Plug in Time')
xlabel('Hour')
ylabel('SOC')
figure(4)
plot(T,PEVT1)
title('Load Profile EV Car')
xlabel('Hour')
ylabel('Load (KW)')
figure(5)
plot(T,PEVT2)
title('Load Profile EV Motorcycle')
xlabel('Hour')
ylabel('Load (KW)')
figure(6)
plot(T,PEVTotal)
title('Load Profile EV')
xlabel('Hour')
ylabel('Load (KW)')