%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input arrival time distribution data
run('\\filestore.soton.ac.uk\users\tsb1e14\mydocuments\MATLAB\Queuing\ReadVecs.m')
Input_ChargePointsVec = [2 3 4 5]; % Number of charging points
Input_MaxGridInput = 150000000000 * 1000; % Maximum power from grid in W
Input_StorageCapacity = 1000 * 1000 * 3600; % Storage capacity of stationary energy storage in J
% Efficiencies
Eff_Grid2EV = 0.95;
Eff_Grid2OVES = 0.95;
Eff_OVES2EV = 0.95;
Eff_OVES = 0.95;
Eff_EV = 0.95;
% Normal distribution of cars using the fast charging station each day
CarsDay_mean = 100;
CarsDay_sd = 10;
% EV charging profile inputs in normal distributions
% EV battery capacity (J)
EV_Bat_cap_mean = 80 * 3600 *1000;
EV_Bat_cap_sd = 8 * 3600 *1000;
% Initial state of charge (%/100)
EV_SOC_init_mean = 0.2;
EV_SOC_init_sd = 0.1;
% Final state of charge (%/100)
EV_SOC_fin_mean = 0.8;
EV_SOC_fin_sd = 0.1;
% SOC when constant power stops (%/100)
EV_SOC_const_mean = 0.25;
EV_SOC_const_sd = 0.025;
% Constant power value (W)
EV_P_const_mean = 400 * 1000;
EV_P_const_sd = 40 * 1000;
% Factor to determine shape of charging curve
EV_k2_mean = 2.0;
EV_k2_sd = 0.1;
certainty = 0.1;
AverageWaitVec = zeros(length(Input_ChargePointsVec),3);
ProbGreat5Vec = zeros(length(Input_ChargePointsVec),3);
iterations_hold = zeros(length(Input_ChargePointsVec),1);
n_loop_count = 100;
for tt = 1:length(Input_ChargePointsVec)
Input_ChargePoints = Input_ChargePointsVec(tt);
for ee = 1:20
wait_temp = zeros(3,1);
grt5_temp = zeros(3,1);
for ff = 1:3
for xx = 1:n_loop_count
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Making inputs and predefining variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Randomly generate the number of EVs using the fast charging station
% during the day
Input_CarsDay = round(normrnd(CarsDay_mean,CarsDay_sd));
% Create a vector of the time of day of arrival of each car
ArrivalTimes = randsample(TimeOfDaySeconds, Input_CarsDay, true, ProbArrivSeconds);
% Create vectors which are the length of the number of cars that use the
% fast charging station during the day
EV_Bat_cap = normrnd(EV_Bat_cap_mean,EV_Bat_cap_sd,Input_CarsDay,1);
EV_SOC_init = normrnd(EV_SOC_init_mean,EV_SOC_init_sd,Input_CarsDay,1);
EV_SOC_init(EV_SOC_init<0.01)=0.01;
EV_SOC_fin = normrnd(EV_SOC_fin_mean,EV_SOC_fin_sd,Input_CarsDay,1);
EV_SOC_fin(EV_SOC_fin>1)=1;
EV_SOC_const = normrnd(EV_SOC_const_mean,EV_SOC_const_sd,Input_CarsDay,1);
EV_P_const = normrnd(EV_P_const_mean,EV_P_const_sd,Input_CarsDay,1);
EV_P_const(EV_P_const>Input_MaxGridInput*Eff_Grid2EV)=Input_MaxGridInput*Eff_Grid2EV;
EV_k2 = normrnd(EV_k2_mean,EV_k2_sd,Input_CarsDay,1);
% Create a matrix, each row is each car and the 100 columns are the maximum
% power that can be delivered to the EV at each state of charge from 1% SoC
% to 100% SoC
EV_max_charge = zeros(Input_CarsDay,100);
TSoC01 = zeros(Input_CarsDay,100);
% Predefine a vector to calculate how long each EV takes to charge
t_charge = zeros(Input_CarsDay,1);
% Calculate the charge power profile of each EV assuming that there are no
% grid constraints, the charge profile is versus state of charge
% Go through each EV
for y = 1:Input_CarsDay
% Go through each state of charge
for x = 1:1000
% If the state of charge is within the constant power section
if x/1000 < EV_SOC_const(y)
% The power is the constant power
EV_max_charge(y,x) = EV_P_const(y);
% Time to change SoC by 0.1%
TSoC01(y,x) = ((0.1 / 100) * EV_Bat_cap(y)) / EV_max_charge(y,x);
% Else if the state of charge is in the decaying section
else
% Calculate the maximum power
EV_max_charge(y,x) = (EV_P_const(y)/exp(-EV_k2(y)*EV_SOC_const(y))) * exp (-EV_k2(y) * x/1000);
% Time to change SoC by 0.1%
TSoC01(y,x) = ((0.1 / 100) * EV_Bat_cap(y)) / EV_max_charge(y,x);
end
end
% Calculate the charging time for each EV
t_charge(y) = sum( TSoC01(y,round(EV_SOC_init(y)*1000): round(EV_SOC_fin(y)*1000)) );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
% Time that constant power is applied for (s)
if EV_SOC_const(y) > EV_SOC_init(y)
t_const = (EV_SOC_const(y)-EV_SOC_init(y))*EV_Bat_cap(y)/(EV_P_const(y)*Eff_EV);
else
t_const = 0;
end
% Constant in curve equation (J)
k3 = EV_Bat_cap(y);
% Constant in curve equation (kW)
k1 = (EV_P_const(y)*Eff_EV)/exp(-EV_k2(y)*EV_SOC_const(y));
% Constant in curve equation
C = ((exp(EV_SOC_const(y)*EV_k2(y)))/(-EV_k2(y)))+(k1*t_const/k3);
% Calculate the length of time charging takes (s)
t_charge(y) = (C-((exp(EV_SOC_fin(y)*EV_k2(y)))/(-EV_k2(y))))*(k3/k1);
%}
end
% State the end time without waiting
ProjectedEnd = ArrivalTimes + t_charge;
StateOfCharge = EV_SOC_init;
% OVES state of charge vectors, at time zero it is 1.0 initially
OVESSoC = 1.0;
OVESSoCVec = zeros(86400,1);
% Create a vector for the cars that are in the queue
QueueVec = [];
CarsCurrentlyCharging = [];
NumCarsCharging = 0;
ChargeStart = NaN(Input_CarsDay,1);
ChargeEnd = NaN(Input_CarsDay,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Go through each second of the day
for i = 1:86400
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Go through all cars in the queue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If there are cars in the queue see if they can start charging
% If there are no cars in the queue do nothing
if isempty(QueueVec)
% Else if there are cars in the queue
else
% Determine how many charging spaces are avaiable
NumSpaces = Input_ChargePoints - NumCarsCharging;
% If there are no spaces do nothing
if NumSpaces == 0;
% Else if there are spaces
else
% Calculate the number of cars that start charging from the
% queue
if NumSpaces >= length(QueueVec)
NumStartCharging = length(QueueVec);
else
NumStartCharging = NumSpaces;
end
% Go through the number of cars that start charging
for j = 1:NumStartCharging
% Set the charging start time of that car to the current
% time
ChargeStart(QueueVec(j)) = i;
% Add the car to cars currently charging
CarsCurrentlyCharging = [CarsCurrentlyCharging QueueVec(j)];
end
% Increase the number of cars charging by the number of cars
% that started charging
NumCarsCharging = NumCarsCharging + NumStartCharging;
% Remove the cars that started charging from the queue
QueueVec = QueueVec(NumStartCharging+1:length(QueueVec));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Go through all cars that arrive in the second
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CarsArriveInSec = find(ArrivalTimes==i);
% If no cars arrive in this second do nothing
if isempty(CarsArriveInSec)
% If cars do arrive then start them charging if possible
else
% Go through all cars that arrive in the second
for j = 1:length(CarsArriveInSec)
% If there is space to start charging
if NumCarsCharging < Input_ChargePoints
% Set the charge start to now
ChargeStart(CarsArriveInSec(j)) = i;
% Add the car to cars currently charging
CarsCurrentlyCharging = [CarsCurrentlyCharging CarsArriveInSec(j)];
% Add one to the number of cars charging
NumCarsCharging = NumCarsCharging+1;
% If there is no space add them to the queue
else
QueueVec = [QueueVec CarsArriveInSec(j)];
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Go through charge added to each car
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If there are no cars are charging all power is available to charge
% the OVES
if NumCarsCharging < 1
PowerOVES = Input_MaxGridInput * Eff_Grid2OVES;
% else for the case if there are cars charging
else
% define a vector for the SoC Added
SoCAdded = zeros(NumCarsCharging,1);
% find the state of charge of the cars currently charging to the
% nearest integer
SoCofCarsCharging = round(1000*(StateOfCharge(CarsCurrentlyCharging)));
MaxChargeCarCurrentCharging = zeros(NumCarsCharging,1);
% find the maximum charge rate of all cars charging
for maxxx = 1:NumCarsCharging
MaxChargeCarCurrentCharging(maxxx) = (EV_max_charge(CarsCurrentlyCharging(maxxx),SoCofCarsCharging(maxxx))) / Eff_EV;
end
% If the required charge rate can be provided just by the grid
if sum(MaxChargeCarCurrentCharging) <= Input_MaxGridInput * Eff_Grid2EV
% the state of charge added to each car is the maximum charge rate
SoCAdded = (MaxChargeCarCurrentCharging * Eff_EV) ./ EV_Bat_cap(CarsCurrentlyCharging);
% The additional power from the grid can b used to charge the OVES
PowerOVES = (Input_MaxGridInput - ((sum(MaxChargeCarCurrentCharging)) / Eff_Grid2EV)) * Eff_Grid2OVES;
% If the power cannot be provided by the grid alone
else
% if the OVES is not empty the power can be provided by the
% grid and the OVES
if OVESSoC > 0
% The cars are therefore charged at the maximum rate
SoCAdded = (MaxChargeCarCurrentCharging * Eff_EV) ./ EV_Bat_cap(CarsCurrentlyCharging);
% The power from the OVES is equal to the maximum rate
% minus the grid connection
PowerOVES = ((Input_MaxGridInput * Eff_Grid2EV) - sum(MaxChargeCarCurrentCharging)) / Eff_OVES2EV;
% else if the OVES is empty
else
% go through each EV and see if there is enough power
% available from the grid to charge it
% initially create a variable for grid power and the number
% charging at full power
PowerAvailFromGrid = Input_MaxGridInput * Eff_Grid2EV;
NumFullPower = 0;
for fullpowercheck = 1:NumCarsCharging
if MaxChargeCarCurrentCharging(fullpowercheck) < PowerAvailFromGrid;
NumFullPower = NumFullPower + 1;
PowerAvailFromGrid = PowerAvailFromGrid - MaxChargeCarCurrentCharging(fullpowercheck);
SoCAdded(fullpowercheck) = ((MaxChargeCarCurrentCharging(fullpowercheck)) * Eff_EV) / EV_Bat_cap(CarsCurrentlyCharging(fullpowercheck));
else
break
end
end
% The remaining grid power is given to the next EV
SoCAdded(NumFullPower+1) = (PowerAvailFromGrid * Eff_EV)/ EV_Bat_cap(CarsCurrentlyCharging(NumFullPower+1));
% There is no power available to give the OVES
PowerOVES = 0;
end
end
% Check to see if any EVs finish charging
% add the state of charge to each EV
StateOfCharge(CarsCurrentlyCharging) = StateOfCharge(CarsCurrentlyCharging) + SoCAdded;
CarFinishCount = 0;
CarsFinish = 0;
% go through each EV charging and see if it finishes charging
for finishCheck = 1:NumCarsCharging
if StateOfCharge(CarsCurrentlyCharging(finishCheck)) > EV_SOC_fin(CarsCurrentlyCharging(finishCheck))
CarsFinish(CarFinishCount+1) = CarsCurrentlyCharging(finishCheck);
CarFinishCount = CarFinishCount + 1;
end
end
% if any cars finish charging
if CarFinishCount > 0
for cenddd = 1:CarFinishCount
ChargeEnd(CarsFinish(cenddd)) = i;
end
% remove cars from the number of cars that finish
NumCarsCharging = NumCarsCharging - CarFinishCount;
% Remove the cars that finish charging from the cars currently
% charging
for removecar = 1: CarFinishCount
CarsCurrentlyCharging = CarsCurrentlyCharging(CarsCurrentlyCharging~=CarsFinish(removecar));
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change the SoC of the OVES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Limit the OVES maximum state of charge to 1
if (OVESSoC + ((PowerOVES * Eff_OVES)/(Input_StorageCapacity))) > 1
PowerOVES = ((1-OVESSoC)*(Input_StorageCapacity)) / Eff_OVES;
end
% Increase the OVES SoC by the power delivered
OVESSoC = OVESSoC + ((PowerOVES * Eff_OVES) /(Input_StorageCapacity));
% Limit the minimum OVES SoC to 0
if OVESSoC < 0
OVESSoC = 0;
end
% Create a vector for the OVES SoC with time
OVESSoCVec(i) = OVESSoC;
end
% If at the final minute any cars are still charging it is assumed that
% they do not have to wait, a warning is displayed if there is a queue at
% the end
for j = 1:length(CarsCurrentlyCharging)
ChargeEnd(CarsCurrentlyCharging(j)) = ChargeStart(CarsCurrentlyCharging(j))+t_charge(CarsCurrentlyCharging(j));
end
if isempty(QueueVec)
else
WarningStillCarsInQueue = 1
end
WaitingTimes = ChargeEnd-ProjectedEnd;
WaitingTimes(WaitingTimes<0)=0;
waitingtime = sum(WaitingTimes);
greatthan5 = length(find(WaitingTimes>300));
%{
Avewaithistcounts = Avewaithistcounts + histcounts(WaitingTimes,bins_for_avewait);
% The time at each C-Rate defined by bins_for_Crate
Time_at_Crate = histcounts(nonzeros(OVESCRateVec),bins_for_Crate);
% Energy through the off-vehicle energy store
Energy_Through_OVES = Input_StorageCapacity * sum(OVESCRateVec(OVESCRateVec>0)) / 3600;
Depths_of_Discharge = zeros(5,1);
% Find the number of times the depth of discharge goes below 20, 40, 60, 80
% 100%
for k = 1:4
Depths_of_Discharge(k)=(length(find(diff((OVESSoCVec-(k*0.2))>0)~=0)+1))/2;
end
Depths_of_Discharge(5)=(length(find(diff((OVESSoCVec-(0.99))>0)~=0)+1))/2;
Depths_of_Discharge = Depths_of_Discharge - [0; Depths_of_Discharge(1:4)];
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wait_temp(ff,1) = wait_temp(ff,1) + waitingtime;
grt5_temp(ff,1) = grt5_temp(ff,1) + greatthan5;
end
end
% largest difference
large_dif = max([abs(wait_temp(1,1)-wait_temp(2,1)),abs(wait_temp(1,1)-wait_temp(3,1)),abs(wait_temp(3,1)-wait_temp(2,1))]);
if large_dif < certainty*min(wait_temp)
break
end
n_loop_count = n_loop_count*2;
clearvars -except grt5 wait Input_CarsDay Input_ChargeTime Input_ChargePoints ProbArrivSeconds TimeOfDaySeconds n_loop Input_MaxGridInput Input_StorageCapacity AverageWaitVec ProbGreat5Vec tt ss n_loop_count wait_temp grt5_temp ee ff Input_ChargePointsVec iterations_hold Input_ChargePointsVec certainty EV_Bat_cap_mean EV_Bat_cap_sd EV_SOC_init_mean EV_SOC_init_sd EV_SOC_fin_mean EV_SOC_fin_sd EV_SOC_const_mean EV_SOC_const_sd EV_P_const_mean EV_P_const_sd EV_k2_mean EV_k2_sd CarsDay_mean CarsDay_sd Eff_Grid2EV Eff_Grid2OVES Eff_OVES2EV Eff_OVES Eff_EV
end
iterations_hold(tt) = n_loop_count;
% wait = mean(wait_temp);
% grt5 = mean(grt5_temp);
AverageWait = wait_temp / (Input_CarsDay*(n_loop_count));
ProbGreat5 = 100*grt5_temp/ (Input_CarsDay*(n_loop_count));
AverageWaitVec(tt,:) = AverageWait;
ProbGreat5Vec(tt,:) = ProbGreat5;
end