Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
Qje=0.36;(*viz kapitola 1.1.1*)
Qjc=8.315;(*viz kapitola 1.1.1*)
Qj=Qje+Qjc;(*viz (8)*)
QTp=QTe+QTc;(*viz (7)*)
cTe=(1/0.92)-(Qje/6);(*viz (18)*)

Vnitřní charakteristika turbosoustrojí QTc = f (Pp) dle dat výrobce

Data={1,14.224,2,19.776,3,25.326,4,30.876,5,36.427,6,41.977}//Flatten;(*data od výrobce turbíny tepelný výkon dodávaná z protitlaku v závislisti na elektrickém výkonu Tab. 7*)
TepelneVykony=Partition[Data,2];
graf=ListPlot[TepelneVykony];
fce=a*(x^b)+8.315;(*nalezení funkce vnitřní charakteristiky*)
find=FindFit[TepelneVykony,fce,{a,b},x];
QTc=fce/.find;(*funkce vnitřní charakteristika*)
graf2=Plot[QTc,{x,0,6},Axes->{True},AspectRatio->1,PlotStyle->{Red,Thickness[0.005]},AxesLabel->{"Pp (MW)","QTc (MW)"},GridLines->Automatic];
Show[graf2,graf];

η = f (Pp) dle výrobce

Data1={1,6.96,2,9.8,3,11.32,4,12.3,5,12.96,6,13.45}//Flatten;
Ucinnost=Partition[Data1,2];
grafU=ListPlot[Ucinnost];
fceU=a*(x^b)+f;
findU=Quiet@FindFit[Ucinnost,fceU,{a,b,f},x];
η=fceU/.findU;
graf2U=Plot[η,{x,0,6},Axes->{True},AspectRatio->1,PlotStyle->{Red,Thickness[0.005]},AxesLabel->{"Pp (MW)","η (%)"},GridLines->Automatic];
Show[graf2U,grafU];

Spotřeba tepla na elektrický výkon QTe=f(Pp) dle dat výrobce

QTe=cTe*x+Qje;(*tepelný výkon odpovídající elektrickému výkonu viz (15)*)
GrafQTe=Plot[QTe,{x,0,6},Axes->{True},AspectRatio->1,PlotStyle->{Red,Thickness[0.005]},AxesLabel->{"Pp (MW)","QTe (MW)"},GridLines->Automatic];

Vnější charakteristika turbosoustrojí QTp=f(Pp) dle dat výrobce

GrafQTp=Plot[QTp,{x,0,6},Axes->{True},AspectRatio->1,PlotStyle->{Red,Thickness[0.005]},AxesLabel->{"Pp (MW)","QTp (MW)"},GridLines->Automatic];

Charakteristiky poměrných přírůstků

cTp=D [QTp,x];(*viz vztah 16*) 
cTc=D [QTc,x];(*viz vztah 17*)

Korekce na skutečné parametry provozovaných soustrojí

ΔcTc=0.1309;(* měrná odchylka od vnitřní měrné spotřeby dané výrobcem(MW/MW) Tab.8*)
ΔcTp=0.15;(* měrná odchylka od vnější měrné spotřeby dané výrobcem(MW/MW)Tab.8*)
ΔcTe=ΔcTp-ΔcTc;(* měrná odchylka od měrné spotřeby tepla na elektrický výkon dané výrobcem(MW/MW)Tab.8*)
Δη=0.133;

Korigované charakteristiky turbosoustrojí

(*Korekce dat naměřených výrobcem na skutečné hodnoty z provozu*)
nQTc=QTc-ΔcTc*x;
nQTp=QTp-ΔcTp*x;
nQTe=nQTp-nQTc;
nη=η+Δη*x;

Vnitřní charakteristika protitlakového trubosoustrojí korigovaná podle skutečných hodnot z provozu

Vnitrni=Plot[nQTc,{x,0,6},Axes->{True},AspectRatio->1,PlotStyle->{Red,Thickness[0.005]},AxesLabel->{"Pp (MW)","QTc (MW)"},GridLines->Automatic,AxesOrigin->{0,0}];

Vnější charakteristika protitlakového trubosoustrojí korigovaná podle skutečných hodnot z provozu;

Vnejsi=Plot[{nQTp},{x,0,6},Axes->{True},AspectRatio->1,PlotStyle->{Red,Thickness[0.005]},AxesLabel->{"Pp (MW)","QTp(MW)"},GridLines->Automatic,AxesOrigin->{0,0}];

Model provozu

ΔIkotelK3=3113.27;(*kJ/kg pro 435°C P3,5 Mpa a H2O 45°C*)
ΔIkotelK45=2588.31;(*kJ/kg pro 0,9 Mpa a H2O 45°C*)
PK45max=5.752; (*maximální celkový výkon kotle K4 a K5 pro ΔIkotelK45 MW*)
QTechnologie=73.56;(* tepelný výkon odebíraný techonologií*)
PsvorkyMin=5; (*elektrický výkon trvale odebíraný technologií*)
nQTc1=nQTc/.x->(Subscript[P, 1]);
nQTc2=nQTc/.x->(Subscript[P, 2]);
nQTp1=nQTp/.x->(Subscript[P, 1]);
nQTp2=nQTp/.x->(Subscript[P, 2]);
nQTc11=2*nQTc1/.(Subscript[P, 1])->(0.5*x);(*vnitřní charakteristika pro paralelně pracující stejně zatížená trubosoustrojí*)
nQTp11=2*nQTp1/.(Subscript[P, 1])->(0.5*x);
(*vnější charakteristika pro paralelně pracující stejně zatížená trubosoustrojí*)
PsvorkyMax=2*FindRoot[nQTc1==0.5*QTechnologie,{(Subscript[P, 1]),5}][[1,2]](*maximální generovaný elektrický výkon při kterém je zajištěn odběr emisního tepla*)
PK3MaxProPelMax=2*nQTp1/.(Subscript[P, 1])->(0.5*PsvorkyMax)(*výkon kotle K3 pro PsvorkyMax MW*)
PK3vykonostPelMax=(PK3MaxProPelMax/ΔIkotelK3)*3600(*výkonost kotle K3 pro PsvorkyMax t/h*)
KotelMin=(QTechnologie-nQTc11+nQTp11)/.x->5;(*minimální výkon kotle K3 pro pokrytí spotřeby tepla technologie a PsvorkyMin MW*)

Navrhovaný výkon kotle v t/h

PnNavrzeno=98;
(*navržený výkon kotle v t/h*)
PnNavrzenoMW=(PnNavrzeno/3600)*ΔIkotelK3;
(*navržený výkon kotle v MW*)

Předem nasmouvaný výkon - optimalizované QTc a QTp pro dosažitelný výkonový rozsah

(*při dodávce nasmlouvané energie se volí rozložení výkonů tak aby byl co nenjnižší tepelný příkon tzn. nejnižší tepelné ztráty*)

RozsahP={{1,6},{1,6}};(*rozsah výkonu turbíny TP 6*)
TeploDodavkaTG={{nQTc1},{nQTc2}};(*minimalizovaná funkce, protitlakový tepelný výkon*)
Psvorky={5,5.5,6,6.005,6.5,7,7.5,8,8.5,9,9.5,10,PsvorkyMax};(*pozadovany vykon v case, od minima technologie po dosažitelné maximum*)


fcePenale[{{a_,b_},i_}]:=Abs[((((Subscript[P, i])-((a+b)/2)))^5)]((1/2)(-Tanh[1000000000((Subscript[P, i])-a)])+0.5+(1/2)Tanh[1000000000((Subscript[P, i])-b)]+0.5)*500000;
fcePenaleProMin[i_]:=fcePenale[{RozsahP[[i]],i}];
fcePenaleProMinSKomb[komb_]:=komb.Table[fcePenaleProMin[i],{i,Length[TeploDodavkaTG]}];

fceNaklady[i_]:=TeploDodavkaTG[[i]]
fceNakladySKomb[komb_]:=komb.Table[fceNaklady[i],{i,1,Length[TeploDodavkaTG]}]
nakladyDiagramSkomb[Pzadany_,komb_]:=(((Evaluate[komb.Table[(Subscript[P, i]),{i,1,Length[TeploDodavkaTG]}]]-Pzadany))^2)*100000000;
fceSoucetproMin[Pzadany_,komb_]:=fcePenaleProMinSKomb[komb]+fceNakladySKomb[komb]+nakladyDiagramSkomb[Pzadany,komb]

pocProMin[komb_]:=Table[komb[[i]]*{(Subscript[P, i]),100},{i,1,Length[TeploDodavkaTG]}]
pocProMinSkomb[komb_]:=Cases[pocProMin[komb],{_Subscript,_}]
fceFindmin[Pzadany_,komb_]:=Module[{vysl},
vysl=FindMinimum[fceSoucetproMin[Pzadany,komb],pocProMinSkomb[komb]];
{vysl}
]
kombAll=Drop[Tuples[{0,1},2],1];
fceHodinaOptimal[hodina_]:=Module[{data},
data=Table[fceFindmin[Psvorky[[hodina]],kombAll[[i]]],{i,1,Length[kombAll]}]//Quiet;
Sort[data][[1]]
](*optimálizované rozdělení výkonů podle požadovaného celkového Pel*)

Průběhy tepelných výkonů v závislosti na elektrickém výkonu

QRedukciS=Table[{Psvorky[[i]], QTechnologie-fceHodinaOptimal[i][[1,1]]},{i,Length[Psvorky]}];(*tepelný výkon předávaný bypassem redukční stanici*)

QTc12=Table[{ Psvorky[[i]],fceHodinaOptimal[i][[1,1]]},{i,Length[Psvorky]}];(*celkový protitlakový tepelný výkon*)

QTp12=Quiet@Table[{Psvorky[[i]],Plus[nQTp1/.fceHodinaOptimal[i][[1,2,1]],If[NumericQ[nQTp2/.fceHodinaOptimal[i][[1,2,2]]],nQTp2/.fceHodinaOptimal[i][[1,2,2]],0]]},{i,Length[Psvorky]}];(*celkový protitlakový tepelný příkon*)

PelVykony1=Table[fceHodinaOptimal[i][[1,2,1,2]],{i,Length[Psvorky]}];(*výkon TG1 v závislosti na požadovaném celekovém Pel*)

PelVykony2=Quiet@Table[If[NumericQ[fceHodinaOptimal[i][[1,2,2,2]]],fceHodinaOptimal[i][[1,2,2,2]],0],{i,Length[Psvorky]}];(*výkon TG2 v závislosti na požadovaném celekovém Pel*)

PK3K4K5=Table[{QRedukciS[[i,1]],QRedukciS[[i,2]]+QTp12[[i,2]]},{i,Length[Psvorky]}];(*celkový výkon dodávaný v páře*)

DodavkaDoSite=ListLinePlot[{QTc12,QRedukciS,QTp12,PK3K4K5},Filling->Bottom,PlotLegends->{"QTc","Výkon redukční stanice","QTp","Výkon praních kotů"},AxesLabel->{"Pp (MW)","Q (MW)"},GridLines->Automatic];
Limit90t=Plot[{If[PnNavrzenoMW>nQTp11,PnNavrzenoMW+PK45max,PnNavrzenoMW]},{x,5,PsvorkyMax},PlotLegends->{"Limit výkonu 90 t/h"},AxesLabel->{"Pp (MW)","Q (MW)"},GridLines->Automatic,PlotStyle->{Black,Thickness[0.007]},AxesOrigin->{5,0}];
DodavkaSlimitem=Show[DodavkaDoSite,Limit90t]

Rozdělení elektricých výkonů TG1 a TG2 v závislosti na požadovaném celkovém výkonu - režim trvalé dodávky

PrubehyPel=ListLinePlot[{Table[{i-1,Drop[Psvorky,{4}][[i]]},{i,Length[Psvorky]-1}],Table[{i-1,Drop[PelVykony1,{4}][[i]]},{i,Length[Psvorky]-1}],Table[{i-1,Drop[PelVykony2,{4}][[i]]},{i,Length[Psvorky]-1}]},Filling->Bottom,PlotLegends->{"TG 1+2","TG 1","TG 2","Pel MAX","Pel MIN"},AxesLabel->{"-","Pp (MW)"},GridLines->Automatic,PlotStyle->{{Green,Thickness[0.007]},{Yellow,Thickness[0.007]},{Blue,Thickness[0.007]}}];
PMINMAX=ListLinePlot[{{{0,PsvorkyMax},{11,PsvorkyMax}},{{0,5},{11,5}}},PlotLegends->{"Pel MAX","Pel MIN"},PlotStyle->{{Red,Thickness[0.007]},{Black,Thickness[0.007]}}];
RozlozeniZatizeni=Show[PrubehyPel,PMINMAX]

Diagram při poskytování podpůrných služeb

VykonyPPS=Plot[{nQTc11,nQTp11,0.5x,QTechnologie-nQTc11,QTechnologie-nQTc11+nQTp11,If[PnNavrzenoMW>nQTp11,PnNavrzenoMW+PK45max,PnNavrzenoMW]},{x,5,PsvorkyMax},PlotLegends->{"QTp","QTc","TG 1 = TG 2","Výkon redukční stanice","Potřebný parní výkon","Navržený výkon + K4 a K5"},AxesLabel->{"Pp (MW)","Q (MW)"},GridLines->Automatic,PlotStyle->{{Black,Thickness[0.007]},{Red,Thickness[0.007]},{Blue,Thickness[0.007]},{Yellow,Thickness[0.007]},{Orange,Thickness[0.007]},{Green,Thickness[0.01]}},AxesOrigin->{5,0}]

ZmenaPpS=Plot[{QTechnologie-nQTc11+nQTp11},{x,5,PsvorkyMax},PlotLegends->{"Výkon parních kotlů pro Pel 5-10,371 MW"},AxesLabel->{"Pp (MW)","Q (MW)"},GridLines->{{5,PsvorkyMax},{KotelMin,PK3MaxProPelMax}},PlotStyle->{{Orange,Thickness[0.007]},{Green,Thickness[0.01]}},AxesOrigin->{5,78}]

ZmenaPkotle=D[{QTechnologie-nQTc11+nQTp11},x];
ZmenaPkotle/.x->6

TypPpS=5;(*RMZt5 - minuty pro uvolnění zálohy*)

SjednanyP=PsvorkyMax-PsvorkyMin;(*velikost regulační zálohy*)

PozadovaRychlostRegulace=ZmenaPkotle*(SjednanyP/TypPpS)/.x->5(*požadovaná rychlost změny výkonu kotle pro poskytnutí snednané zálohy ve sjednaném čase*)

Pprocentne=(PozadovaRychlostRegulace/PnNavrzenoMW)*100
(*maximální růst výkonu v % jmenovitého výkonu*)

KontrolaP[v_]:=If[(PnNavrzenoMW>nQTp11/.x->v),"Výkon kotle je dostatečný pro požadovaný elektrický výkon"v"MW","Výkon kotle je nedostatečný pro požadovaný elektrický výkon"v"MW"](*kontrola výkonu kotle*)

KontrolaP[9.5]

Export["zatizeni.png",ZmenaPpS];