%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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