% Modelování energetických toků ES
clc;
close all;


% % Simulace bytového ES
% load('simulace_1.mat');
% 
% byt = zeros(24, 86400);
% dum = zeros(3,86400);
% fve = vyroba(33.59,3);
% for i = 1:24
%     byt(i,:) = model_spotreby_1;
%     dum(mod(i,3)+1,:) = dum(mod(i,3)+1,:) + byt(i,:);
% end
% 
% metrika_S3(byt(1,:),'Byt č.1', 15);
% % 
% metrika_S3(dum - fve,'RDS',15);
% metrika_S1(dum - fve,'RDS',15);

% % Simulace různorodého ES
% 
% % load('simulace_2.mat') 
% % 
% oblacnost = randn(1,48*3).*randn(1,48*3);
% 
% RD1 = model_spotreby_1(4.5,3)- vyroba(5.25,3,200,40,oblacnost);
% RD2 = model_spotreby_1(4,3);
% RD3 = model_spotreby_1(2.5,3) - vyroba(11.81,3,105,55,oblacnost);
% RD4 = model_spotreby_1(5.2,3)- [vyroba(1.31,1,180,0,oblacnost); zeros(2,86400)];
% RD5 = model_spotreby_1(3.5,3)- [vyroba(2.62,2,240,30,oblacnost); zeros(1,86400)];
% 
% metrika_S3(RD1,'RD1');
% metrika_S3(RD2,'RD2');
% metrika_S3(RD3,'RD3');
% metrika_S3(RD4,'RD4');
% metrika_S3(RD5,'RD5');
% metrika_S3(RD1+RD2+RD3+RD4+RD5,'RDS');
% metrika_S1(RD1,'RD1');
% metrika_S1(RD2,'RD2');
% metrika_S1(RD3,'RD3');
% metrika_S1(RD4,'RD4');
% metrika_S1(RD5,'RD5');
% metrika_S1(RD1+RD2+RD3+RD4+RD5,'RDS');
   


function energie = metrika_S3(vykon_s, nazev, agr_min, faze)  % Sekundový výkon v W, časové okno v min (default 15), výstup v kWh, selektivní zpracování dat
arguments
    vykon_s {mustBeNumeric};
    nazev {mustBeText} = 'A';
    agr_min = 15;
    faze {mustBeInteger} = 1;
end
if mod(length(vykon_s),agr_min*60) ~= 0
    error("Délka měřicího okna není soudělná s délkou měřeného úseku.");
end
[f, ~] = size(vykon_s);
if f == 1
    vykon_s(mod(faze+2,3)+1,:) = vykon_s;
    vykon_s(mod(faze,3)+1,:) = 0;
    vykon_s(mod(faze+1,3)+1,:) = 0;
end
mw_s = agr_min*60;
energie = zeros(2,length(vykon_s)/mw_s);
    for i = 1:length(vykon_s)
        for ii = 1:3
            if vykon_s(ii,i)>0
                energie(1,ceil(i/mw_s)) = energie(1,ceil(i/mw_s)) + vykon_s(ii,i);
            else
                energie(2,ceil(i/mw_s)) = energie(2,ceil(i/mw_s)) + vykon_s(ii,i);
            end
        end
    end
stair1=zeros(2,length(vykon_s));
for i=1:length(vykon_s)
    stair1(1,i) = energie(1,ceil(i/mw_s))/(1000*mw_s);
    stair1(2,i) = energie(2,ceil(i/mw_s))/(1000*mw_s);
end
fig1 = figure;
fig1.Name = 'Metrika S3';
fig1.Units = 'centimeters';
fig1.Position=[2 2 19 10];
fig1.Resize = 'off';
en = stairs(0:1/3600:86399/3600, [vykon_s(:,:)'/1000 stair1(:,:)']);
en(1).Color = '#666';
en(2).Color = '#999';
en(3).Color = '#ccc';
en(4).Color = '#f00';
en(5).Color = '#090';
en(4).LineWidth = 1;
en(5).LineWidth = 1;
grid on
title(append('Průběh změřený měřicím zařízením \Sigma3 v bodě ', nazev, ', agregační interval ', int2str(agr_min),' min'));
ylabel("P_p (kW)");
xlabel("t (h)");
legend('L1', 'L2','L3',"Odběr", "Dodávka", 'location', 'southeast')
xlim([0 24]);
xticks([0:2:24]);
% 
% fprintf("\Sigma3 metrika, měřicí okno %.0f min: Energie spotřebovaná/dodaná v bodě \n %s: %.3f/%.3f kWh \n",mw_min, nazev, sum(energie(1,:))/(3600000), sum(energie(2,:))/(3600000));
end

function energie = metrika_S1(vykon_s, nazev, agr_min, faze)  % Sekundový výkon v W, časové okno v min (default 15), výstup v kWh, součtové zpracování dat
arguments
    vykon_s {mustBeNumeric};
    nazev {mustBeText} = 'B';
    agr_min = 15;
    faze {mustBeInteger} = 1;
end
[f, ~] = size(vykon_s);
if f == 1
    vykon_s(mod(faze+2,3)+1,:) = vykon_s;
    vykon_s(mod(faze,3)+1,:) = 0;
    vykon_s(mod(faze+1,3)+1,:) = 0;
end
mw_s = agr_min*60;
energie = zeros(2,length(vykon_s)/mw_s);
    for i = 1:length(vykon_s)
        if sum(vykon_s(1:3,i))>0
            energie(1,ceil(i/mw_s)) = energie(1,ceil(i/mw_s)) + sum(vykon_s(1:3,i));
        else
            energie(2,ceil(i/mw_s)) = energie(2,ceil(i/mw_s)) + sum(vykon_s(1:3,i));
        end
    end
% stair1=zeros(2,length(vykon_s));
% for i=1:length(vykon_s)
%     stair1(1,i) = energie(1,ceil(i/mw_s))/(1000*mw_s);
%     stair1(2,i) = energie(2,ceil(i/mw_s))/(1000*mw_s);
% end
% fig1 = figure;
% fig1.Name = 'Metrika S1';
% fig1.Units = 'centimeters';
% fig1.Position=[2 2 19 10];
% fig1.Resize = 'off';
% en = stairs(0:1/3600:86399/3600, [vykon_s(:,:)'/1000 stair1(:,:)']);
% en(1).Color = '#666';
% en(2).Color = '#999';
% en(3).Color = '#ccc';
% en(4).Color = '#f00';
% en(5).Color = '#090';
% en(4).LineWidth = 1;
% en(5).LineWidth = 1;
% grid on
% title(append('Průběh změřený měřicím zařízením \Sigma1 v bodě ', nazev, ', agregační interval ', int2str(agr_min),' min'));
% ylabel("P_p (kW)");
% xlabel("t (h)");
% legend('L1', 'L2','L3',"Odběr", "Dodávka", 'location', 'southeast')
% xlim([0 24]);
% xticks([0:2:24]);
% 
% fprintf("\Sigma1 metrika, měřicí okno %.0f min: Energie spotřebovaná/dodaná v bodě \n %s: %.3f/%.3f kWh \n",mw_min, nazev, sum(energie(1,:))/(3600000), sum(energie(2,:))/(3600000));
end

function vykon = m2_to_kW(plocha)                       % funkce převádějící plochu na dostupný výkon fve
vykon=0.16*1.64*floor(plocha/1.64);
end

function Spotreba=model_spotreby_1(Spotreba_rok, pocet_fazi)        % Vstup MWh/rok, -, výstup ve W. Stěžejní funkce, která vrátí denní diagram odběru zadaného množství domácností po 1 min. Vykreslí diagram spotřeby a hodnoty skutečného a očekávaného hodinového průměru.
arguments
    Spotreba_rok {mustBeNumeric} = 1.5;                 % MWh/rok, byt ~1.5, RD ~4.5
    pocet_fazi {mustBeInteger} = 1;
end

%% zavedení parametrů modelu
tdd = [0.42, 0.40,0.38,0.37,0.40,0.46,0.50,0.53,0.55,0.56,0.57,0.56,0.56,0.57,0.58,0.62,0.70,0.76,0.77,0.75,0.68,0.58,0.48,0.44,0.42];
t = 1/3600:1/3600:24;    
obalka = spline(0:1:24,tdd, t);
Pocet_zapnuti_den= [500,  20,     70,     2];          % Průměrný počet zapnutí za den
Pp_min =           [2,    90,     400,    600];           % Minimální příkon spotřebiče v dané kategorii (kW)
Pp_max =           [25,   200,    1000,   1000];          % Maximální příkon spotřebiče v dané kategorii (kW)
T_zap_min =        [1,    1,      0.2,    20];             % Minimální délka cyklu (min)
T_zap_max =        [60,   5,      5,      45];           % Maximální délka cyklu (min)

for i = 1:pocet_fazi
    Spotreba(i,:) = 0*t+50+90*round(rand-0.15*i);                                     % vektor spotreby (W)
end

for i = 1:4
    Pocet_zapnuti = (Pocet_zapnuti_den(i) - random('Exponential',1))*Spotreba_rok/8.7;
    f = 0;
    while true
        f = f+1;
        if Pocet_zapnuti <= 0
            break
        end
        Pp_sp = probexp(Pp_min(i), Pp_max(i));
        T_zap = probexp(T_zap_min(i)*60, T_zap_max(i)*60);
        bias = 0;
%         if i == 3 || i == 4
%             bias = Pp_sp/Pp_max(i);
%         end
        t_start = probability_function(obalka,bias)+length(obalka);
        for ii = 1:T_zap
            if Pp_sp > 250
                Pp_sp_on = Pp_sp * (0.9*rand()+0.1);
            else
                Pp_sp_on = Pp_sp;
            end 
            Spotreba(mod(f,pocet_fazi)+1,mod(t_start-round(T_zap*1/3)+ii,length(Spotreba))+1) = Spotreba(mod(f,pocet_fazi)+1,mod(t_start-round(T_zap*1/3)+ii,length(Spotreba))+1) + Pp_sp_on;
        end
        Pocet_zapnuti = Pocet_zapnuti - random('Exponential',1);
    end
end  
% fprintf('Spotreba: %f MWh/rok',sum(Spotreba*365/3600000000,'all'));
end

function Spotreba = model_spotreby_2(Spotreba_rok, pocet_fazi) % Model 2, vrací denní diagram odběru o zadané průměrné spotřebě
arguments
    Spotreba_rok {mustBeNumeric} = 1.5;                 % MWh/rok, byt ~1.5, RD ~4.5
    pocet_fazi {mustBeInteger} = 1;
end
tdd = [0.42, 0.40,0.38,0.37,0.40,0.46,0.50,0.53,0.55,0.56,0.57,0.56,0.56,0.57,0.58,0.62,0.70,0.76,0.77,0.75,0.68,0.58,0.48,0.44,0.42];
t = 1/3600:1/3600:24;    
obalka = spline(0:1:24,tdd, t);

base_load = Spotreba_rok/0.09 * (randn(1,1)/6+1);     % W, byt (10% průměrné denní spotřeby ~ 4 kWh/den)
Spotreba_den = Spotreba_rok*1000 * (randn(1,1)/6+1);
Prikon = [5 2200];                              % Rozmezí příkonu
Doba_zapnuti = [60 60*440];                     % Rozmezí doby zapnutí
Spotreba = zeros(1,length(t))+base_load;
Energie_za_den = base_load*24;

while 1
    Prikon_sp = random('Beta',0.4,4)*(Prikon(2)-Prikon(1))+Prikon(1);
    Doba_zap_sp = round((Doba_zapnuti(2)-Doba_zapnuti(1))*(5+25/2200)*(1./Prikon_sp-1/2200)+Doba_zapnuti(1));
    if abs(Energie_za_den + Prikon_sp*Doba_zap_sp/3600 - Spotreba_den) > abs(Energie_za_den - Spotreba_den)
        break 
    end
    Cas_zap_sp = probability_function(obalka,0.2);
    for ii = 1:Doba_zap_sp
        Spotreba((mod(round(Cas_zap_sp+ii-Doba_zap_sp/3+length(Spotreba)),length(Spotreba))+1)) = Spotreba((mod(round(Cas_zap_sp+ii-Doba_zap_sp/3+length(Spotreba)),length(Spotreba))+1)) + Prikon_sp;
    end
    Energie_za_den = Energie_za_den + Prikon_sp*Doba_zap_sp/3600;
end
end

function vyroba = vyroba(vykon, pocet_fazi, azimut, sklon, oblacnost)      % Stěžejní funkce vracející časový průběh výroby (W) s náhodnou oblačností pro pocet výroben se zadanými výkony v kW a azimuty 90(východ) - 270(západ)
arguments
    vykon (1,:) {mustBeFloat}
    pocet_fazi (1,1) {mustBeInteger} = 3
    azimut (1,:) {mustBeFloat} = 180;
    sklon (1,:) {mustBeFloat} = 40;
    oblacnost (1,:) {mustBeNumeric} = randn(1,48*3).*randn(1,48*3);
end
t_s = 1/3600:1/3600:24;

for i = 2:length(oblacnost)
    oblacnost(i) = oblacnost(i) + oblacnost(i-1);
end
oblacnost = spline(1:48*3, oblacnost, 0:48*3/length(t_s):48*3-1/length(t_s));
oblacnost = oblacnost/6;
oblacnost = oblacnost*0.8/mean(oblacnost);
for i = 1:length(oblacnost)
    if oblacnost(i) > 1
        oblacnost(i) = 1;
    elseif oblacnost(i) < 0.1
        oblacnost(i) = 0.1;
    end
end

doba_svitu = 14;                     % h, doba od nuly do nuly výkonu
osvit = cos((t_s-12)*pi/(doba_svitu)).^2;
for i = 1:length(t_s)
    if t_s(i) < (24-doba_svitu)/2 || t_s(i) > 12+doba_svitu/2
        osvit(i) = 0;
    end
end
osvit_vyrobna = zeros(1,length(t_s));
osvit_vyrobna = cos((sklon-35)/60)*vykon*cos((t_s-((sklon/90)*azimut/15+12*(1-sklon/90)))*pi/(24-sklon/7.5));
osvit_vyrobna = osvit_vyrobna.*oblacnost;
for ii = 1:length(t_s)
    if osvit_vyrobna(ii) < 0.1*vykon
        osvit_vyrobna(ii) = 0.1*vykon;
    end
end
osvit_vyrobna = osvit_vyrobna.*osvit;
if azimut < 90 || azimut > 270
    warning("Pozor: orientace některé FV výrobny přesahuje interval 90° - 270°, výsledky mohou být zkreslené. \n");
end

vyroba = zeros(pocet_fazi, length(osvit_vyrobna));
for i = 1:pocet_fazi
    vyroba(i,:) = 1000*osvit_vyrobna/pocet_fazi;
end
end

%% Pomocné funkce

function out = probexp(low, top)                        % Interní funkce vracející hodnotu mezi minimem a maximem podle exponenciálního rozdělení
    out = random('Exponential',top);
    while low > out || out > top
        out = random('Exponential',top);
    end
end

function output = probability_function(vect, bias)      % Interní funkce vracející pozici ve vstupním vektoru vect s pravděpodobností danou hodnotami v tomto vektoru. Předpokládá kladné hodnoty. 
arguments
    vect (1,:) {mustBeNumeric}
    bias (1,1) {mustBeNumeric} = 0;
end
Min = min(vect);
vect=vect-Min*bias;
s = sum(vect);
    if s == 0
        error("Vektorem musí být kladná nenulová funkce!")
    end
    roll = rand()*s;
    i = 1;
    while roll > 0
        roll = roll - vect(i);
        i = i +1;
    end
    output = i;
end