Quiet@Remove["Global`*"];
$HistoryLength=2;
Off[FindRoot::lstol]
Off[NIntegrate::izero]
SetDirectory[NotebookDirectory[]];

c=4186.;(*měrná tepelná kapacita vody*)
(*známý stav pro podlahové vytápění*)
znamyStav={Ptopny-> 6000,T2->T1-5,T1->45, Tin->20 };
rceZnamyStav={m*c*(T1-T2)==Ptopny,Ptopny==k*(((T1+T2)/2)-Tin)}//.znamyStav;
(*vypočtené konstatní parametry*)
konstParametry=Union[Solve[rceZnamyStav][[1]],{Tin->20} ]

(*dosazení uložených závislostí*)
{Ptopny[Tvzduch_],PtopnyR[Tvzduch_],PtopnyC[Tvzduch_]}=<<"Ptop";
{PelHP[Tvzduch_,Tvoda_],PtopHP[Tvzduch_,Tvoda_],COP[Tvzduch_,Tvoda_]}=<<"NIBE AMS10-8";

(*Rovnice systému kdy čerpadlo pracuje na hranici bivalence*)
rceHranice[Tvzduch_]={m*c*(T1-T2)==Ptopny[Tvzduch],Ptopny[Tvzduch]==k*(((T1+T2)/2)-Tin),COP[Tvzduch,T1]*Pel==Ptopny[Tvzduch],Ptopny[Tvzduch]==PtopHP[Tvzduch,T1]}/.konstParametry;
(*počáteční hodnoty od kterých FindRoot se hledá řešení*)
poc={{Pel,2000},{T1,40},{T2,40},{Tvzduch,-5}};
(*řešení rovnic a nalezení Teploty bivalence*)
resHranice=FindRoot[rceHranice[Tvzduch],poc]
TmezBivalence=Tvzduch/.resHranice;
Print["Teplota bivalence je ",TmezBivalence," °C"]

(*Rovnice systému kdy čerpadlo dokáže pokrýt 100% tepelných ztrát*)
rce100[Tvzduch_]={m*c*(T1-T2)==Ptopny[Tvzduch],Ptopny[Tvzduch]==k*(((T1+T2)/2)-Tin),COP[Tvzduch,T1]*Pel==Ptopny[Tvzduch]}/.konstParametry;
(*řešení rovnic*)
res100[Tvzduch_]:=Union[FindRoot[rce100[Tvzduch],{{Pel,2000},{T1,40},{T2,40}}],{PtopOstatni->0}];

(*Rovnice systému kdy čerpadlo pracuje v bivalenci*)
rceBivalentni[Tvzduch_]={m*c*(T1-T2)==PtopHP[Tvzduch,T1],PtopHP[Tvzduch,T1]==k*(((T1+T2)/2)-Tin),COP[Tvzduch,T1]*Pel==PtopHP[Tvzduch,T1],PtopOstatni+PtopHP[Tvzduch,T1]==Ptopny[Tvzduch]}/.konstParametry;
Cases[rceBivalentni[-5],_Symbol,{0,∞}]//Union;
Length/@{%,rceBivalentni[-5]};
resBival[Tvzduch_]:=FindRoot[rceBivalentni[Tvzduch],{{Pel,2000},{T1,40},{T2,40},{PtopOstatni,1000}}];

(*sjednocení řešení systému s tepelným čerpadlem*)
res[Tvzduch_]:=Quiet@res100[Tvzduch]/;Tvzduch≥TmezBivalence;
res[Tvzduch_]:=Quiet@resBival[Tvzduch]/;Tvzduch<TmezBivalence;

(*nastavení grafů*)
ps={AxesLabel->{"t[s]","T[°C]"},PlotLabel->"  Průměrná denní teplota jako funkce času",PlotStyle->{Thickness[0.007],RGBColor[1,0,0]},GridLines->Automatic};

(*elektrický příkon tepelného čerpadla*)
Pelektricky[Tvzduch_]:=Pel/.res[Tvzduch];
datPel={#,Pelektricky[#]}&/@Range[-15,20,0.2];
Pelint=Interpolation[datPel]; (*interpolovaný*)
(*topný výkon druhého zdroje pro bivalenci k TČ*)
u=0.899; (*účinnost zdroje*)
Post[Tvzduch_]:=(PtopOstatni/u)/.res[Tvzduch];
datPost={#,Post[#]}&/@Range[-15,20,0.2];
Postint=Interpolation[datPost];(*interpolovaný*)
(*celkový topný výkon druhého zdroje*)
Postcint[Tvzduch_]=Postint[Tvzduch]+((PtopnyR[Tvzduch])/u);(*interpolovaný*)

Plot[{Pelint[T],Postcint[T]},{T,-15,20},AxesLabel->{"Tv [°C]","P [W]"},PlotStyle->{{Thickness[0.007],RGBColor[0,0,1]},{Thickness[0.007],RGBColor[1,0,0]}}, PlotLegends->{"Elektrický příkon TČ","Celkový příkon plynového kotle"},GridLines->{Union[Table[i,{i,-20,25,5}],{TmezBivalence}],Automatic}]

(*dostupná data průměrné měsíční teploty pro danou lokalitu*)
teplotaData={{1,-2.5},{2,-1.2},{3,2.4},{4,6.7},{5,11.7},{6,15},{7,16.5},{8,15.8},{9,12.5},{10,7.8},{11,2.4},{12,-1},{13,-2.5}};
pomfce[month_]=Interpolation[teplotaData][month];
prumRocniTeplota=7.2;
rok=365.*24*3600; (*rok v sekundách*)
denrok=365 ;(*rok ve dnech*)
(*vytvoření smyčky*)
teplotaMesicni[m_]:=pomfce[m]/;1≤m<13;
teplotaMesicni[m_]:=teplotaMesicni[m-12]/;m≥13;
teplotaMesicni[m_]:=teplotaMesicni[m+12]/;m<1;
teplotaSec[t_]:=teplotaMesicni[(t/(rok/12))];
teplotaDen[t_]:=teplotaMesicni[(t/(denrok/12))];
prumTeplotaGraf =(1/denrok)NIntegrate[teplotaDen[t],{t,0,denrok}];
prumTeplota=(1/rok)NIntegrate[teplotaSec[t],{t,0,rok}];
(*Závislost teploty na čase v sekundách*)
teplota[t_]:=teplotaSec[t]-prumTeplota+prumRocniTeplota;
(*Závislost teploty na čase ve dnech*)
teplotaGraf[t_]:=teplotaDen[t]-prumTeplotaGraf+prumRocniTeplota;
(*teplota  při které je zahájena a ukončena topná sezóna*)
touttopim=13; 

Plot[{teplota[t],touttopim},{t,0,2*rok},AxesLabel->{"t[s]","T[°C]"},PlotStyle->{{Thickness[0.005],GrayLevel[0]},{Thickness[0.007],RGBColor[1,0,0]}},GridLines->Automatic]
castopenipocatek=t/.FindRoot[teplota[t]==touttopim,{t,2*(10^7),3*(10^7)}];
castopenikonec=t/.FindRoot[teplota[t]==touttopim,{t,4*(10^7),5*(10^7)}];

pocatektopnesezony =Floor[(castopenipocatek/(24*3600))];
konectopnesezony=Floor[(castopenikonec/(24*3600))]-365;

Plot[{teplotaGraf[t],touttopim},{t,0,365},AxesLabel->{"t[den v roce]","T[°C]"},PlotLabel->"Průměrná teplota jako funkce času",PlotStyle->{{Thickness[0.005],GrayLevel[0]},{Thickness[0.007],RGBColor[1,0,0]}},GridLines->{Union[Table[i,{i,0,365,30}],{konectopnesezony,pocatektopnesezony}],Table[i,{i,-10,25,2.5}]},Ticks->{Union[Table[i,{i,0,365,30}],{konectopnesezony,pocatektopnesezony}],Union[Table[i,{i,-10,25,2.5}],{touttopim}]},PlotRange->{{0,365},{-5,20}},AspectRatio->1/1.5,ImageSize->Scaled[0.8]]

Print["Počátek topné sezóny =  ",pocatektopnesezony,"  den v roce"]
Print["Konec topné sezóny =  ",konectopnesezony,"  den v roce"];

(*výpočet potřeby tepla*)
potrebaTepla=NIntegrate[Ptopny[teplota[t]],{t,castopenipocatek,castopenikonec}];(*okruh s podl. vytápěním*)
potrebaTeplaR=NIntegrate[PtopnyR[teplota[t]],{t,castopenipocatek,castopenikonec}];(*okruh s radiatory*)
potrebaTeplaC=NIntegrate[PtopnyC[teplota[t]],{t,castopenipocatek,castopenikonec}];(*celkem*)

(*výpočet spotřeby energie*)
spotrebaElektriny=NIntegrate[Pelint[teplota[t]],{t,castopenipocatek,castopenikonec}];(*tepelného čerpadla*)
spotrebaOstatni=NIntegrate[Postint[teplota[t]],
{t,castopenipocatek,castopenikonec}];(*druhého zdroje pro bivalenci*)
spotrebaOstatniR =potrebaTeplaR/u;(*druhého zdroje pro okruh s radiátory*)
spotrebaOstatniC=NIntegrate[Postcint[teplota[t]],{t,castopenipocatek,castopenikonec}];(*druhého zdroje celkem*)
spotrebaEnergie =spotrebaOstatniC+spotrebaElektriny;(*celková spotřeba energie*)

(*výpočet topného faktoru*)
prumernyCOPc=(potrebaTeplaC/spotrebaEnergie);
prumernyCOP=(potrebaTepla/(spotrebaElektriny+spotrebaOstatni));

(*data pro výpis výsledků*)
dataGJ ={potrebaTepla,potrebaTeplaR,potrebaTeplaC,spotrebaElektriny,spotrebaOstatni,spotrebaOstatniR,spotrebaOstatniC,spotrebaEnergie}*(1/(10^9));
datakWh=dataGJ *((10^9)/(3600*1000));
dataText = {"pro okruh s podl. vytápěním " ,"pro okruh s radiátory ","CELKEM ","tepelného čerpadla ","druhého zdroje pro bivalenci ","druhého zdroje pro okruh s radiátory ","druhého zdroje celkem ","CELKEM "};

(*vypsání a uložení výsledků*)
For[i=1,i<=Length[dataGJ],i++,
If[i==1,Print["---POTŘEBA TEPLA---"]];
If[i==4,Print["---SPOTŘEBA ENERGIE---"]];
Print[dataText[[i]],dataGJ[[i]]," GJ, tj. ",datakWh[[i]]," kWh"];
]
Print["---TOPNÝ FAKTOR---"];
Print["průměrný COP TČ je ",prumernyCOP," [-]"];
Print["účinnost druhého zdroje je ",u*100," [%]"];
Print["průměrný COP vytápění je ",prumernyCOPc," [-]"];
Print["---SPOTŘEBA ENERGIE BEZ TČ---"];
Print[dataGJ[[3]]/u, " GJ, tj. ",datakWh[[3]]/u," kWh"];
Put[{spotrebaElektriny,spotrebaOstatniC,potrebaTeplaC/u},"energy"]