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