Off[General::"spell1"];

Tin = 25;
TextNom = -12;
kW = 1000. ;
PzNom = 10  kW;
reskDum = Solve[
            PzNom == 
              kDum*(Tin - 
                    TextNom)][[1]];
c2 = 4186;
resPom = Solve[{PzNom == 
                  kTopeni*(((T1nom + T2nom)/2) - Tin), 
                m2*c2*(T1nom - T2nom) == PzNom} /. {T1nom -> 
                  55, T2nom -> 
                  50}][[1]];
dosStav = Union[reskDum, resPom];

c1 = 1000;
m1 = 0.5;
rceCerpadlo = {Qin == m1*c1*(Tv - Tout), 
              Qout == m2*c2*(Tk - T2), Pel + Qin == Qout, 
              Qout == ϵ*Pel} /. ϵ -> 
              0.5 (( Tk + 273)/(Tk - Tv)) /. 
          dosStav;
rceKotel = {(T1 - Tk)*m2*c2 == 
              Pk, (T1 - T2)*m2*c2 == 
              kTopeni*(((T1 + T2)/2) - Tin), 
            kTopeni*(((T1 + T2)/2) - Tin) == kDum*(Tin - Texter)} /. 
          Union[{Pk -> 0}, dosStav];
Δ1 = Tzem - Texter;
Δ2 = Tzem - Tv;
kS = 150;
Tzem = 2;
rceZemniVymenik = {Qvymenik == m1*c1*(Tv - Texter), 
          Qvymenik == 
            kS*((Δ1 + Δ2)/2)};

ClearAll[rceDohr];
rceDohr[Texter_] = 
      Union[rceKotel, rceCerpadlo, 
        rceZemniVymenik]
Length[rceDohr[5]]
Union[
      Cases[rceDohr[
          5], _Symbol, {0, ∞}]]
% // Length

res[
          Teplotavenku_] := Solve[
            rceDohr[Teplotavenku]][[1]];
prikon[Teplotavenku_] := 
        Pel /. res[Teplotavenku];
pl1 = 
        Plot[prikon[tt]/1000, {tt, -12, 13}, 
          PlotStyle -> {Thickness[0.008]}, 
          AxesLabel -> {"Tvenku[°C]", "P[kW]"}];

PmaxCerpadlo = 8000;

TmezBivalence = 
        Texter /. Solve[
              rceDohr[Texter] /. 
                Qout -> 
                  PmaxCerpadlo][[1]];
Print["hranice bivalence je ", 
      TmezBivalence, " °C"]

rceKotelSBivalenci = {(T1 - Tk)*m2*c2 == 
              Pk, (T1 - T2)*m2*c2 == 
              kTopeni*(((T1 + T2)/2) - Tin), 
            kTopeni*(((T1 + T2)/2) - Tin) == kDum*(Tin - Texter)} /. 
          Union[dosStav];

rceCerpadloSBivalenci = {Qin == strida*m1*c1*(Tv - Tout), 
              Qout == strida*m2*c2*(Tk - T2), 
              Pel + Qin == Qout, 
              Qout == ϵ*Pel} /. ϵ -> 
              0.5 (( Tk + 273)/(Tk - Tv)) /. dosStav;

rceDohrSBivalenci[5] = 
      Union[rceKotelSBivalenci, rceCerpadloSBivalenci, 
        rceZemniVymenik]

rceVysl[
          Texter_] := (rceDohrSBivalenci[Texter] /. {Pk -> 0, 
                strida -> 1}) /; 
          Texter ≥ TmezBivalence;
rceVysl[
          Texter_] := (rceDohrSBivalenci[Texter] /. {strida -> 1, 
                Qout -> PmaxCerpadlo}) /; 
          Texter < TmezBivalence;

rceVysl[Texter_] := 
        Union[Solve[(rceDohrSBivalenci[Texter] /. {Pk -> 0, 
                      strida -> 1})][[1]], {Pk -> 0}] /; 
          Texter ≥ TmezBivalence;
rceVysl[Texter_] := 
        Union[Solve[(rceDohrSBivalenci[Texter] /. {strida -> 1, 
                      Qout -> 
                        PmaxCerpadlo})][[1]], {Qout -> 
                PmaxCerpadlo}] /; Texter < TmezBivalence;

prikonSBivalenci[Texter_] := 
        Pk + Pel /. rceVysl[Texter];
vykonSBivalenci[Texter_] := 
        Pk + Qout /. rceVysl[Texter];
pl2 = 
        Plot[prikonSBivalenci[tt]/1000, {tt, -12, 13}, 
          PlotStyle -> {Thickness[0.008], Hue[0]}, 
          AxesLabel -> {"Tvenku[°C]", "P[kW]"}];
Show[pl1, pl2]

\(\(pl2 = 
        Plot[vykonSBivalenci[tt]\/prikonSBivalenci[tt], {tt, \(-12\), 13}, 
          PlotStyle -> {Thickness[0.008], Hue[0]}, 
          AxesLabel -> {"Tvenku[°C]", "ϵCelk[-]"}];\)\)

TemperatureData = {{1, -1.4}, {2, - .2}, {3, 3.7}, {4, 
            8.9}, {5, 13.8}, {6, 17.2}, {7, 19}, {8, 18.5}, {9, 14.8}, {10, 
            9.1}, {11, 
            4.1}, {12,  .3}, {13, -1.4}};
proklad = 
        Interpolation[TemperatureData, 
          InterpolationOrder -> 3];
ClearAll[tepVenku];
tepVenku[mesic_] := 
        proklad[mesic] /; 
          1 ≤ mesic < 13;
tepVenku[mesic_] := 
        tepVenku[mesic - 12] /; 
          mesic ≥ 13;
tepVenku[mesic_] := tepVenku[mesic + 12] /; mesic < 1;

Off[NIntegrate::"slwcon"];
rok = 365.25*24*3600;
mesicSekund = (1/12)*rok;
teplota[tsekund_] := 
        tepVenku[(tsekund/mesicSekund)];
prumernaTeplota = 
        (1/rok)*NIntegrate[
            teplota[t], {t, 0, rok}];
zadanaPrumernaTeplota = 6;
teplVenDef[t_] := 
        teplota[t] - prumernaTeplota + zadanaPrumernaTeplota;

Plot[{teplVenDef[t], 13}, {t, 0.8, 1.8  rok}, 
        PlotStyle -> {{Hue[0], Thickness[0.008]}, {Hue[0.6], 
              Thickness[0.008]}}];

tpoc = 
      t /. FindRoot[
          teplVenDef[t] == 13, {t, 2.5*(10^7)}]
tkon = 
      t /. FindRoot[
          teplVenDef[t] == 13, {t, 4.5*(10^7)}]