Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
den=3600*24.;

{{Pel,Ptop},info}=<<"struct";
(*Ta vstupni vzduch do TC, Tw vystupni voda z TC*)
Plot3D[{Pel[Ta,Tw],Ptop[Ta,Tw]},{Ta,-40,40},{Tw,20,60}];

(*msteckouMax - hmotnostni prutok mezi TC a aku nadrzi, kVykon reprezentuje velikost TC, Maccu - hmotnost vody v nádrži v kg, Taccu teplota ukum. nadrze, msteckou prutok v kg/s*)
(*funkce vydejDers vraci derivaci teploty v aku nadrzi a derivaci odebrane elekricke energie*)

ClearAll[vydejDers];
MAglobal=8000.;
vydejDers[{Ta_,Taccu_}]:=Module[{Maccu=MAglobal,c=4186.,ΔBezp=3,msteckouMax=0.6,Twmin,msteckou,Twpom,Tw,msteckouMinTw,derTaccu,derEel,kVykon=0.99},
Twmin=Pel[[1,2,1]]+ΔBezp;
Twpom=Tw/.Quiet@FindRoot[msteckouMax*c*(Tw-Taccu)==kVykon*1000*Ptop[Ta,Tw],{Tw,20}];
msteckouMinTw=kVykon*((1000*Ptop[Ta,Twmin])/(c*(Twmin-Taccu)));
derTaccu:=kVykon*((1000*Ptop[Ta,Twpom])/(c*Maccu))/;Twpom≥Twmin;
derTaccu:=kVykon*((1000*Ptop[Ta,Twmin])/(c*Maccu))/;Twpom<Twmin;
derEel:=kVykon*1000*Pel[Ta,Twpom]/;Twpom≥Twmin;
derEel:=kVykon*1000*Pel[Ta,Twmin]/;Twpom<Twmin;
{derTaccu,derEel}
];

datDers=Quiet@Flatten[Table[{Ta,Taccu,vydejDers[{Ta,Taccu}]},{Ta,-30,40,N[Pi/E]},{Taccu,4,60,N[E/Pi]}],1];

intDers=Interpolation[datDers];
{intDers[-5,10],vydejDers[{-5,10}]}

posun=1*365.25*24*3600;

teploty=(<<"teplotyBrno");

Take[teploty,8];

teplDenniPom=Interpolation[teploty];
tmaxPom=teplDenniPom[[1,1,2]];
(*resamplovani teplot*)
teplotypom={#,teplDenniPom[#]}&/@Range[0.,tmaxPom,1800];
teplDenni1=Interpolation[teplotypom];
(*tady je zohlednen posun v casove souradnici*)
teplDenni=teplDenni1[#+posun]&;

secondsInMonth=365.*24*3600/12;
Plot[teplDenni[t],{t,0,12*secondsInMonth}];

ClearAll[vydej2];
vydej2[{{TStart_,Tzadana_},{tStart_,Tvzduch_}}]:=Module[{c=4186.,MA=MAglobal,derivace,RungeKuttaKrok,Δt=10.*60,krok,T1a,sol,TAa,T1aa,ta,tb,tEnd,Etopna,ErealEl,COPReal,TAb,Eelb,Eela},

derivace[{t_,stav:{TAccu_,Eel_}}]:=intDers[Tvzduch[t],TAccu];
RungeKuttaKrok[{t_,y_}]:=Module[
{k1,k2,k3,k4},
k1=derivace[{t,y}];
k2=derivace[{t+Δt/2,y+(Δt/2)*k1}];
k3=derivace[{t+Δt/2,y+(Δt/2)*k2}];
k4=derivace[{t+Δt,y+Δt*k3}];
{t+Δt,y+(Δt/6)*(k1+2k2+2k3+k4)}
];
sol=NestWhileList[RungeKuttaKrok,{tStart,{TStart,0}},#[[2,1]]≤Tzadana&];
{{ta,{TAa,Eela}},{tb,{TAb,Eelb}}}=Take[sol,-2];
tEnd=(-ta TAb+TAa tb+ta Tzadana-tb Tzadana)/(TAa-TAb);
Etopna=c*MA*(Tzadana-TStart);
ErealEl=(Eelb TAa-Eela TAb+Eela Tzadana-Eelb Tzadana)/(TAa-TAb);
COPReal=Etopna/ErealEl;
{ErealEl,Etopna,COPReal,{tStart,tEnd}}
];
vydej2[{{10,50},{180000,teplDenni}}]

(*nepripravujeme stale stejne mnozstvi teple vody, ale topime vodu pro budouci spotrebu tepla na vytapeni*)
(*uvazujeme nesetrvacnou budovu s tepelnou ztratou*)
Pzpom[Ta_]=Fit[{{-12,25000},{13,0}},{1,Ta},Ta];
Pz[Ta_]:=Pzpom[Ta]/;Ta≤13;
Pz[Ta_]:=0/;Ta>13;
Plot[Pz[Ta],{Ta,-20,20}];
(*minimalni potrebna teplota vody vstupujici do otopne soustavy aby dokazala soustava predat vykon*)
ClearAll[TwNaVytopeni];
TwNaVytopenipom[Ta_]=Fit[{{-18,55},{13,20}},{1,Ta},Ta];
TwNaVytopeni[Ta_]:=TwNaVytopenipom[Ta]/;Ta≤13;
TwNaVytopeni[Ta_]:=20/;Ta>13;

ClearAll[vydejParametryZitrka];
vydejParametryZitrka[tStart_]:=Module[{nCasu=51,casy,hodnoty,Δt,den=24*3600.,energiePotrebna,prumernyVykon,maxVykon,tEnd,teploty,teplotaPrumerna},
Δt=den/(nCasu-1);
tEnd=tStart+den;
casy=Range[tStart,tEnd,Δt];
teploty=teplDenni/@casy;
hodnoty=Pz[teplDenni[#]]&/@casy;
energiePotrebna=(Δt/3.)*Total[{1.,4.,1}.#&/@Partition[hodnoty,3,2]];
teplotaPrumerna=(Δt/(3.*den))*Total[{1.,4.,1}.#&/@Partition[teploty,3,2]];
prumernyVykon=energiePotrebna/den;
maxVykon=Max[hodnoty];
{tEnd,{If[teplotaPrumerna<13,energiePotrebna,10000],maxVykon,prumernyVykon,teplotaPrumerna}}

];
vydejParametryZitrka[250*den]

(*pro dany start dneska vypocte parametry zitrka a z nich mezi jakymi teplotami vody v nadrzi bychom se meli pohybovat*)
ClearAll[vydejTStartTzadana];
vydejTStartTzadana[tStart_]:=Module[{energiePotrebna,maxVykon,prumernyVykon,tEnd,teplotaPrumerna,TStart,Tzadana,c=4186.},
{tEnd,{energiePotrebna,maxVykon,prumernyVykon,teplotaPrumerna}}=vydejParametryZitrka[tStart];
TStart=TwNaVytopeni[teplotaPrumerna];
Tzadana=TStart+energiePotrebna/(c*MAglobal);
{TStart,Tzadana}

];
vydejTStartTzadana[1*den]

(*funkce energie vypocte potrebnou energii na naukumulovani tepla pro budouci den v case tStart*)

ClearAll[energie];
myR:=Round[100#]/100.&;
energie[tStart_]:=First@Quiet@vydej2[{vydejTStartTzadana[tStart],{tStart,teplDenni}}];

ClearAll[dejBestStart];
dejBestStart[zacatek_]:=Module[{data,dost,minEel,maxEel,dostMax,intData,tmin,tmax,pom,tminEl,tmaxEl,datapom},data={#,energie[#]}&/@(zacatek+Range[0,den,den/24.]);
datapom={#,energie[#]}&/@(zacatek+Range[0,den,den/(4*24.)]);
pom=SortBy[datapom,Last];
{tminEl,minEel}=First[pom];
{tmaxEl,maxEel}=Last[pom];
dost={t->tminEl};
dostMax={t->tmaxEl};
{minEel,data,dost,maxEel,dostMax}
];
dejBestStart[100*den];

(*vrati graf, kde jsou spotreby elektricke energie v optimalnim (modre), pesimalnim (cervene) a pevne zadanem (zelene) case v hodinach*)
dejGrafy[zacatek_,hodinaFix_Integer]:=Module[{minEel,data,dost,maxEel,dostMax,Efix},
{minEel,data,dost,maxEel,dostMax}=dejBestStart[zacatek];
Efix=data[[hodinaFix+1,2]];
Show[ListPlot[data/.{t_,Ee_}:>{(t-zacatek)/3600,Ee/(10^6)},AxesLabel->{"tStart [hour]","Eel [MJ]"},PlotLabel->"Emin= "<>ToString[myR@(minEel/(10^6))]<>" MJ"<>"   Emax= "<>ToString[myR@(maxEel/(10^6))]<>" MJ"<>"   E"<>ToString[hodinaFix]<>"= "<>ToString[myR@(Efix/(10^6))]<>" MJ",PlotStyle->{Black,PointSize[0.01]},AxesOrigin->{0,0},GridLines->Automatic],ListPlot[{{(t-zacatek/.dost)/3600,minEel/(10^6)}},PlotStyle->{Blue,PointSize[0.02]}],ListPlot[{{(t-zacatek/.dostMax)/3600,maxEel/(10^6)}},PlotStyle->{Red,PointSize[0.02]}],ListPlot[{{hodinaFix,Efix/(10^6)}},PlotStyle->{Green,PointSize[0.02]}]]
];
uk11=dejGrafy[60*den,7]
Export["uk11.jpg",uk11]

dejGrafyTepl[zacatek_,hodinaFix_Integer]:=Module[{minEel,data,dost,maxEel,dostMax,Efix,plTeplota,t1,t2},
{minEel,data,dost,maxEel,dostMax}=dejBestStart[zacatek];
Efix=data[[hodinaFix+1,2]];
{t1,t2}={zacatek,zacatek+den};
plTeplota=Plot[teplDenni[t*3600+t1],{t,0,(t2-t1)/3600},PlotStyle->{Magenta,Thickness[0.0044]}];
Show[ListPlot[data/.{t_,Ee_}:>{(t-zacatek)/3600,Ee/(10^9)},AxesLabel->{"tStart [hour]","Eel [GJ]  Ta [°C]"},PlotLabel->"Emin= "<>ToString[myR@(minEel/(10^9))]<>" GJ"<>"   Emax= "<>ToString[myR@(maxEel/(10^9))]<>" GJ"<>"  Epevný "<>ToString[hodinaFix]<>"= "<>ToString[myR@(Efix/(10^9))]<>" GJ",PlotStyle->{Black,PointSize[0.01]},AxesOrigin->{0,0},GridLines->Automatic],ListPlot[{{(t-zacatek/.dost)/3600,minEel/(10^9)}},PlotStyle->{Blue,PointSize[0.02]}],ListPlot[{{(t-zacatek/.dostMax)/3600,maxEel/(10^9)}},PlotStyle->{Red,PointSize[0.02]}],ListPlot[{{hodinaFix,Efix/(10^9)}},PlotStyle->{Green,PointSize[0.02]}],plTeplota,PlotRange->All]

];
uk12=dejGrafyTepl[60*den,7]
Export["uk12.jpg",uk12]

ClearAll[energieFix];
energieFix[indexDne_,hodinaFix_]:=energie[indexDne*den+(hodinaFix)*3600];
myR@(energieFix[103,8]/(10^6))

energieFixRok[hodinaFix_]:=Module[{dataFix},
dataFix=energieFix[#,hodinaFix]&/@Range[1,364];
{hodinaFix,Plus@@dataFix}
];
dataEFix={#,energieFixRok[#]}&/@Range[0,23];
dataEFix>>"dataEFix"

dataDejBestStart={#,dejBestStart[#*den]}&/@Range[1,364];
dataDejBestStart>>"dataDejBestStart";

dataEFix=<<"dataEFix";
dataDejBestStart=<<"dataDejBestStart";

datafix=Last/@dataEFix;
ListPlot[datafix]
intfix=Interpolation[datafix];
minFix=FindMinimum[intfix[h],{h,8}]

maxFix=FindMaximum[{intfix[h],0<h<23},{h,21}];
minBod={h/.minFix[[2]],((minFix[[1]])/(10^9))};
maxBod={h/.maxFix[[2]],((maxFix[[1]])/(10^9))};

pla=ListPlot[datafix/.{h_,E_}:>{h,(E/(10^9))},AxesLabel->{"tStart [hour]","Eel [GJ]"},PlotStyle->{Black,PointSize[0.01]},AxesOrigin->{0,0},GridLines->Automatic,PlotLabel->"EminFix="<>ToString[myR@minBod[[2]]]<>"GJ  EmaxFix="<>ToString[myR@maxBod[[2]]]<>"GJ"];
Show[pla,ListPlot[{minBod},PlotStyle->{Blue,PointSize[0.02]}],ListPlot[{maxBod},PlotStyle->{Red,PointSize[0.02]}]]

minBod

Eoptim={#[[1]],#[[2,1]]}&/@dataDejBestStart;
Length/@{Eoptim,dataDejBestStart};
uk2=plOptim=ListPlot[Eoptim/.{day_,Eel_}:>{day,Eel/(10^6)},AxesLabel->{"day","Emin [MJ]"},GridLines->Automatic,PlotStyle->{Blue,PointSize[0.006]}]
Export["uk2.pdf",uk2]
Export["uk2.png",uk2,ImageResolution->500]

dataPesim;

test[{d_,Eel_}]:=True;
dataPesim={#[[1]],#[[2,-2]]}&/@dataDejBestStart;
Epesim=Select[dataPesim,test];
plPesim=ListPlot[Epesim/.{day_,Eel_}:>{day,Eel/(10^6)},AxesLabel->{"day","Emax [MJ]"},GridLines->Automatic,PlotStyle->{Red,PointSize[0.006]}];
Length/@{Epesim,dataDejBestStart}

Show[plOptim,plPesim];

dataBestFix=Drop[{#,energieFix[#,h/.minFix[[2]]]}&/@Range[1,366],{180}];
plBestFix=ListPlot[dataBestFix/.{day_,Eel_}:>{day,Eel/(10^6)},AxesLabel->{"day","Emax [MJ]"},GridLines->Automatic,PlotStyle->{Green,PointSize[0.006]}];
Length[dataBestFix]

uk3=plBestFixAll=Show[plOptim,plPesim,plBestFix,PlotRange->All]
Export["uk3.pdf",uk3]
Export["uk3.png",uk3,ImageResolution->500]

dataWorstFix={#,energieFix[#,h/.maxFix[[2]]]}&/@Range[1,366];
plWorstFix=ListPlot[dataWorstFix/.{day_,Eel_}:>{day,Eel/(10^6)},AxesLabel->{"day","Emax [MJ]"},GridLines->Automatic,PlotStyle->{Green,PointSize[0.006]}];

plWorstFixAll=Show[plOptim,plPesim,plWorstFix]

EWorstFix=myR@(Plus@@(Last/@dataWorstFix)/(10^9));
EBestFix=myR@(Plus@@(Last/@dataBestFix)/(10^9));
EoptimVys=myR@(Plus@@(Last/@Eoptim)/(10^9));
EpesimVys=myR@(Plus@@(Last/@Epesim)/(10^9));
datVys={EoptimVys,EBestFix,EWorstFix,EpesimVys}
((datVys-EoptimVys)/EoptimVys)*100

minBod

{10.28057938487944`,27.207050187908376`}