ClearAll["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
SetAttributes[exp,{HoldAll,Listable}];
exp[co_]:=Module[{a,graf},a=ToString[Unevaluated[co]];graf:=Evaluate[co];Export[Evaluate[a<>".png"],graf,ImageResolution->300]];

kCVypZap=0.5;
Off[FindMinimum::"fmmp"];
dataPminPmax={{100,625},{80,500},{75,400},{30,300}};
dataAcka={{500,8,0.004},{400,6.4,0.0048},{600,7.9,0.005},{400,7.5,0.0055}};
dataCinCout=kCVypZap*{{3000,1500},{3000,1500},{2000,1000},{1000,800}};
n=Length[dataAcka];
vykony=(Subscript[P, #])&/@(rn=Range[n]);
f[x_]:=(x^2);
pen[x_,a_,b_,δa_,δb_]:=(10^3)Module[{res,vyr,xx,k,ka,kb},k=(10^4)Min[Abs[{(1/a),(1/b)}]];vyr[xx_]:=(ka*f[(xx-a)]+kb*f[(xx-b)])*(2+Tanh[(xx-b)*k]-Tanh[(xx-a)*k]);res=Solve[{vyr[a*(1-0.01δa)]==100,vyr[b*(1+0.01δb)]==100},{ka,kb}][[1]];vyr[x]/.res];
ukpen[x_]=pen[x,20,30,5,5];
plPenale=Plot[ukpen[x],{x,0,50},PlotRange->All,PlotStyle->{Red,Thick},AxesLabel->{"P [MW]","penale"},GridLines->Automatic];
exp[plPenale];

penale[i_]:=pen[vykony[[i]],dataPminPmax[[i,1]],dataPminPmax[[i,2]],5,5];
Plot[penale[1],{(Subscript[P, 1]),0,750},PlotRange->All,PlotStyle->{Red,Thick},AxesLabel->{"P [MW]","penale"},GridLines->Automatic];

cena[i_]:=Module[{Pi},
Pi=vykony[[i]];
{1,Pi,(Pi^2)}.dataAcka[[i]]
];

penaleVykony=Plus@@(penale/@rn);
cenaVari=Plus@@(cena/@rn);
penaleDodavka[Pdodavka_]:=200((((Plus@@vykony)-Pdodavka))^2);
proFindMinimum[Pdodavka_]:=penaleVykony+penaleDodavka[Pdodavka]+cenaVari;

obor[i_]:={vykony[[i]],Mean[dataPminPmax[[i]]]//N};
obory=obor/@rn;

min[Pdodavka_]:=FindMinimum[proFindMinimum[Pdodavka],obory];
cenaHodina[Pdodavka_]:=cenaVari/.min[Pdodavka][[2]];
datazatizeni={600,800,1100,1400,1000,1200,1000,1400,1600,1600,1400,800,400,200,1400,1200,450,300,400};
pom1=min/@datazatizeni//Quiet;
nakladyDen=cenaVari/.(#[[2]]&/@pom1);
Print["celkem naklady= ",Plus@@nakladyDen," -" ]

(*optimalizace denniho provozu pres cely den*)

nk=(2^n)-1;
komb=((Take[#,-4]&)/@(IntegerDigits[Range[nk],2]/.{q__Integer}:>{0,0,0,q}))//Union;

nakladyHodina[Pzadany_,kombinace_]:=Module[{penMinMax,penDodavka,kpd=100,kpmn=10,obors,cen,promin,min,m,ps},
penMinMax=kpmn*(penale/@rn).kombinace;
penDodavka=kpd*((((vykony.kombinace)-Pzadany))^2);
obors=Cases[Part[obory,kombinace*rn],{_,_}];
cen=(cena/@rn).kombinace;
promin=cen+penDodavka+penMinMax;
min=FindMinimum[promin,obors]//Quiet;
{cen/.(m=min[[2]]),kombinace,m,ps={penMinMax,penDodavka}/.m,Plus@@ps}
];

vydejHodina[Pzadany_]:=Module[{dats},
dats=nakladyHodina[Pzadany,#]&/@komb;
Take[#,3]&/@Sort@Select[dats,#[[-1]]<1&]
];

cenaprepnuti[v1_,v2_]:=Sum[If[Sign[(v1-v2)][[i]]<0,dataCinCout[[i,1]],If[Sign[(v1-v2)][[i]]>0,dataCinCout[[i,2]],0]],{i,n}];

start={vydejHodina[datazatizeni[[1]]][[1]],1};

krok[{{c_,komb_,dosP_},n_}]:=Module[{zatizeni,mozne,sprepnutim},
zatizeni=datazatizeni[[n+1]];
mozne=vydejHodina[zatizeni];
sprepnutim=Sort[mozne/.{ci_Real,kombi_List,dosPi_}:>{ci+1*cenaprepnuti[komb,kombi],kombi,dosPi}][[1]];
{sprepnutim,n+1}
];

dats=NestList[krok,start,Length[datazatizeni]-1];
cenaDenSVypZap=Plus@@(#[[1,1]]&/@dats);
Print["celkem denni naklady= ",cenaDenSVypZap," -"]

pom=dats/.{{ci_Real,kombi_List,dosPi_},nh_}:>{nh,ci,kombi,dosPi,datazatizeni[[nh]],(ci/(datazatizeni[[nh]]))};
p1=BarChart[#[[2]]&/@pom,AxesLabel->{"hodina","c [Kc]"}]
p2=BarChart[#[[5]]&/@pom,AxesLabel->{"hodina","P [MW]"}]
p3=BarChart[#[[-1]]&/@pom,AxesLabel->{"hodina","ch [(Kc/MWh)]"}]
exp[{p1,p2,p3}]

proexp=pom/.(Subscript[P, a_]):>"P"<>ToString[a]/.(a_->b_):>ToString[a]<>"="<>ToString[b];
Export["vysledky.xls",proexp]