%Definice vstupních proměnných (technické parametry) BESS
Crate = 1;         % Hodnota vyb. a nab. Crate (NMC, NCA)       [-]
Pmax = 5;          % Maximální výkon/výše poskytnuté rezervy    [MW]
n_system = 0.95;   % Celková účinnost systému (BESS + trafo)    [-]
pom_CP = 1.6;      % Poměr mezi Cin a Pmax (30 min krit.)       [-]
Cin = Pmax*pom_CP; % Nejmenší potřebná jmen. kapacita BESS      [MWh]
SoC_start=0.5*Cin; % Nastavení poč. hodnoty SoC 50%             [MWh]
fmin = -10;        % Nastavení dolní hranice pásma necitlivosti [mHz]
fmax = 10;         % Nastavení horní hranice pásma necitlivosti [mHz]
t_off = 120;       % Čas odpojení od FCP při krit. stavu        [min]
%Načtení historických dat frekvence v ES
f_input_2013 = xlsread ('2013');
f_input_2014 = xlsread ('2014');
f_input_2015 = xlsread ('2015');
f_input_2016 = xlsread ('2016');
f_input_2017 = xlsread ('2017');
%f_input_2018 = xlsread ('2018');
%f_max = xlsread ('day_max_2018');
%f_kor=xlsread ('KOR');
%Model ročního provozu BESS pro FCP
f_input=f_input_2017; % volba historických dat pro model
%f_input=f_kor;
%f_input=f_max;
f_deviance = (f_input-50)*1000;
f_DB=f_deviance;
time=length(f_input); % počet minut v roce
SoC=zeros(time,1);
P_DB=zeros(time,1);
P_DF=zeros(time,1);
P_kor=zeros(time,1);
E_bat=zeros(time,1);
E_fcp=zeros(time,1);
E_kor=zeros(time,1);
P=zeros(time,1);
Pcelk=zeros(time,1);
SoC(1)=SoC_start;
n_kor=0;
n = 0;
while (n < time)
    n = n + 1;
    if (SoC(n)< 0.9*Cin) && (SoC(n)>0.1*Cin) 
       if (f_deviance(n)<= fmax) && (f_deviance(n)>=fmin) 
           if  SoC(n)<= 0.4*Cin && (f_deviance(n)>=0)
               P_DB(n)= -0.1*Pmax*((0.3*Cin)/SoC(n));
           elseif SoC(n)>=0.6*Cin && (f_deviance(n)<=0)
               P_DB(n)= 0.1*Pmax*(SoC(n)/(0.6*Cin));
           else
               P_DB(n)=0;
           end
           Pcelk(n)=P_DB(n);
       else
            if f_DB(n) > 0 
                P(n) = (Pmax/200) * f_DB(n)*(-1);
                if SoC(n)<=0.4*Cin
                    P_DF(n)=(P(n)*1.2)-P(n);
                else
                    P_DF(n)=0;
                end
                Pcelk(n)=P(n)+P_DF(n);
            elseif f_DB(n) < 0
                P(n) = (Pmax/200) * f_DB(n)*-1;
                if SoC(n)>=0.6*Cin
                    P_DF(n)=(P(n)*1.2)-P(n);
                else
                    P_DF(n)=0;
                end
                Pcelk(n)=P(n)+P_DF(n);
            else
                P(n) = 0;
                Pcelk(n)=0;
            end    
       end
       E_fcp(n) = ((P(n)*-1)/60);
       E_bat(n) = ((Pcelk(n)*-1)/60);
       SoC(n+1) = SoC(n)+E_bat(n);
    else 
        if SoC(n)>= 0.9*Cin
            while (SoC(n)>= 0.6*Cin)
            P_kor(n)=(0.3*Cin)/(t_off/60);
            Pcelk(n)=P_kor(n);
            E_kor(n) = ((P_kor(n)*-1)/60);
            E_bat(n) = ((P_kor(n)*-1)/60);
            SoC(n+1) = SoC(n)+E_bat(n);
            n=n+1;
            end
        n=n-1;
        n_kor=n_kor+1;
        elseif SoC(n)<= 0.1*Cin
            while (SoC(n)<= 0.4*Cin)
            P_kor(n)=(0.3*Cin)/(t_off/60)*-1;
            Pcelk(n)=P_kor(n);
            E_kor (n)= ((P_kor(n)*-1)/60);
            E_bat(n) = ((P_kor(n)*-1)/60);
            SoC(n+1) = SoC(n)+E_bat(n);
            n=n+1;
            end
        n=n-1;
        n_kor=n_kor+1;
        end
    end
end
SoC_proc=((SoC)/Cin)*100;
%Aplikace účinnosti systémy
nbat=1/n_system;
%Analýza výsledků ročního provozu BESS pro FCP
%Průměrný SoC během roku
SoC_avg = mean (SoC_proc);
figure
histogram(SoC_proc,20);
%Dodaná energie do sítě v rámci FCP
E_posfcp=sum(E_fcp(E_fcp>0));
E_negfcp=sum(E_fcp(E_fcp<0));
E_fcp=E_posfcp+(E_negfcp*-1);
%Celková energie do sítě v rámci FCP s využitím stupňů volnosti
E_pos=sum(E_bat(E_bat>0));
E_neg=sum(E_bat(E_bat<0));
E=E_pos+(E_neg*-1);
%Energie využita pro balancování SoC
E_pos_soc=E_pos-E_posfcp;
E_neg_soc=E_neg-E_negfcp;
Esoc=E_pos_soc+(E_neg_soc*-1);
%Celkové ztráty energie během provozu BESS během roku
Ez_pos=(E_pos*nbat)-E_pos;
Ez_neg=(E_neg*nbat*-1)+E_neg;
%Energie využitá pro balancování SoC
E_soc=E-E_fcp;
%Celková nakoupená korektivní energie po dosažení lim. hodnot SoC
E_poskor=sum(E_kor(E_kor>0))*nbat;
E_negkor=sum(E_kor(E_kor<0))*nbat;
%Počet stavů dosažení limitních hodnot SoC
n_k=((abs(E_negkor)+E_poskor))/n_kor;
n_korneg=E_negkor/n_k*-1;
n_korpos=E_poskor/n_k;
%Výkon během FCP
figure
histogram(P,50);
%Celkové rozložení výkonu během roku
figure
histogram(Pcelk,40);
%Výpočet jednoho EFC a 80% DoD
EFC=Cin;
DoD=0.8*EFC;
EFC_year_fcp=E/EFC;
EFC_year=E_fcp/EFC;
DoD_year=E/DoD;
figure
plot (f_deviance,Pcelk,'g*');
xlim ([-200 200]);
ylim ([-Pmax Pmax]);
xlabel ('odchylka f [mHz]');
ylabel ('Celkový výkon BESS [MW]');
figure
plot(Pcelk)
xlim ([0 time]);
xlabel ('čas [min]');
ylabel ('Výkon BESS [MWh]');
figure
plot(P_kor)
xlim ([0 time]);
xlabel ('čas [min]');
ylabel ('Korekční výkon [MWh]');
% SoC v procentech
figure
plot(SoC_proc)
xlim ([0 time]);
xlabel ('čas [min]');
ylabel ('SoC BESS [%]');