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

(*v souboru dataPID.csv jsou data ve tvaru {cas,vstupPIDregulatoru=in,vystupPIDregulatoru=out}, predpokladame out=p*in[t]+d*in'[t]+i*(∫_0^t)in[τ]\[DifferentialD]τ. Najdete konstanty p,i,d a napiste (slovne) jak byste overovali svoji metodu. Pokud to overeni naprogramujete, mate bonusove body*)
Quiet@Remove["Global`*"];
f[t_]:=10.*Sin[t]+5*Sin[3t]+3*Sin[5t];
tmax=10;
Plot[f[t],{t,0,tmax}];
p=10;
i=7;
d=0.5;
sol=NDSolve[{int'[t]==f[t],int[0]==0},int,{t,0,tmax}][[1]];
iSlozka[t_]:=i*int[t]/.sol;
dSlozka[t_]=d*f'[t];
pSlozka[t_]=p*f[t];
Plot[dSlozka[t],{t,0,tmax}];
vystup[t_]:=pSlozka[t]+iSlozka[t]+dSlozka[t];
data={#,f[#],vystup[#]}&/@Range[0.,tmax,(tmax/500)];
Export["dataPID.csv",data];

data=Import["dataPID.csv"];
vystup=Interpolation[data/.{t_,in_,out_}:>{t,out}];
vstup=Interpolation[data/.{t_,in_,out_}:>{t,in}];
{tmin,tmax}=vystup[[1,1]];
casy=data/.{t_,in_,out_}:>t;
Δt=casy[[2]]-casy[[1]];
hodnotyIntegralu=Δt*Prepend[Accumulate[data/.{t_,in_,out_}:>in],0]//Drop[#,-1]&;
intIntegral=Interpolation[{casy,hodnotyIntegralu}\[Transpose]];
datproFit={vstup[#],intIntegral[#],vstup'[#],vystup[#]}&/@Range[tmin,tmax,((tmax-tmin)/100)];
Fit[datproFit,{pp,ii,dd},{pp,ii,dd}]

(*V souboru "dataDobeh.csv" jsou data ve tvaru {t,ω[t]} soustroji motor-ventilator. Moment ventilatoru je -(10^(-6))*((ω[t])^2). Urcete moment setrvacnosti soustroji, konstantu treni, na jakou uhlovou rychlost se soustroji rozebehne pri hnacim momentu 1 N*m a za jakou dobu dosahne uhlova rychlost poloviny ustalene hodnoty pri tomto momentu. Popiste, jak byste overovali presnost sveho algoritmu pro zasumela data. Pokud to naprogramujete, extra body.*)
Remove["Global`*"];
J=0.3;
B=0.012;
kv=(10^(-6));
ω0=0.45*314;
tmax=2*60;
sol=NDSolve[{J*ω'[t]==-B*ω[t]-kv*((ω[t])^2),ω[0]==ω0},ω,{t,0,tmax}][[1]];
ω=ω/.sol;
Plot[ω[t],{t,0,tmax},PlotRange->All,AxesOrigin->{0,0}];
data={#,ω[#]}&/@Range[0.,tmax,(tmax/100)];
Export["dataDobeh.csv",data];

data=Import["dataDobeh.csv"];
intω=Interpolation[data];
{tmin,tmax}=intω[[1,1]];
datproFit={intω[#],intω'[#],-kv*((intω[#])^2)}&/@Range[0.,tmax,(tmax/500)];
fit=Fit[datproFit,{omega,deromega},{omega,deromega}]

Bmer=Coefficient[fit,omega];
Jmer=Coefficient[fit,deromega];
Mrozbeh=1;
tmaxRozbeh=200;
sol=NDSolve[{ωr'[t]*Jmer==-Bmer*ωr[t]-kv*((ωr[t])^2)+Mrozbeh,ωr[0]==0},ωr,{t,0,tmaxRozbeh}][[1]];
ωrr[t_]:= ωr[t]/.sol;
Plot[ωrr[t],{t,0,tmaxRozbeh},PlotRange->All]

ωinf=ωrr[tmaxRozbeh];
ωPul=0.5ωinf;
tPul=t/.FindRoot[ωrr[t]==ωPul,{t,0.01tmax}];
Print["J= ",Jmer,"  B= ",Bmer,"  ωinf= ",ωinf,"  tPul= ",tPul]

(*Mate zadane INV, diskont, zivotnost a rocni usporu. Urcete NPV=-INV+(∑_(i=1)^zivotnost)(rocniUspora/(((1+0.01*diskont))^i)), IRR takove, ze 0==-INV+(∑_(i=1)^zivotnost)(rocniUspora/(((1+0.01*irr))^i)) a nejak se popasujte s realnou dobou navratnosti definovanou ta, ze 0==-INV+(∑_(i=1)^Treal)(rocniUspora/(((1+0.01*diskont))^i))*)
Remove["Global`*"];
INV=50.5*(10^6);
diskont=3.4;
rocniUspora=4.1*(10^6);zivotnost=20;
NPV=-INV+(∑_(i=1)^zivotnost)(rocniUspora/(((1+0.01*diskont))^i));

IRR=irr/.FindRoot[0==-INV+(∑_(i=1)^zivotnost)(rocniUspora/(((1+0.01*irr))^i)),{irr,diskont}];

datinv=Table[Reverse@{T,-INV+(∑_(i=1)^T)(rocniUspora/(((1+0.01*diskont))^i))},{T,0,zivotnost}];
Treal=Interpolation[datinv,InterpolationOrder->1][0];
Print["NPV= ",NPV,"  IRR= ",IRR,"  Treal= ",Treal]

(*Dvema vedenimi o R1=10Ω a R2=15Ω se z DC zdroje o U=750kV ma prenest celkem 1000 MW. Urcete minimalni ztraty na vedeni a vykony a napeti na vystupu vedeni*)
Remove["Global`*"];
U=750*(10^3);
R1=10;
R2=15;
Pcelk=1000*(10^6);
u1=U-R1*i1;
u2=U-R2*i2;
P1=u1*i1;
P2=u2*i2;
min=FindMinimum[{R1*(i1^2)+R2*(i2^2),Pcelk==P1+P2},{i1,10},{i2,10}]
{((min[[1]])/Pcelk),{P1,P2}/.min[[2]],{R1*(i1^2),R2*(i2^2)}/.min[[2]],{u1,u2}/.min[[2]]}

lagrange=R1*(i1^2)+R2*(i2^2)+λ*(Pcelk-(P1+P2));
0.==D[lagrange,#]&/@{i1,i2,λ};
sol=Solve[%][[2]]
{u1,u2}/.sol

(*je dana DR {y''(t)+2 y'(t)+y(t)==((sin(t))/((Sqrt[(t^(3/2))+1])+1)),y(0)==0,y'(0)==a}, sestrojte graf zavislosti hodnoty y[t=1] na hodnote a pro a ∈<0,20>*)
(*najdete hodnotu a pro kterou plati y[t=1]=2*)
Remove["Global`*"];
rce[a_]:={y''[t]+2y'[t]+y[t]==((Sin[t])/(1+(Sqrt[(t^(3/2))+1]))),y[0]==0,y'[0]==a};
tmax=1;
amax=20;
sol[a_]:=y[tmax]/.NDSolve[rce[a],y,{t,0,tmax}];
data=Table[{a,sol[a]},{a,0,amax,0.05amax}];
int=Interpolation[data];
Plot[int[a],{a,0,amax}]
Interpolation[Reverse/@data][2]

(*soubor matSoustavy.csv obsahuje rozsirenou matici soustavy, posledni sloupec je vektor pravych stran*)
(*urcete soucet prvku vektoru reseni a cislo podminenosti soustavy s vyuzitim spektralni normy (te co vyuziva vlastni cisla, v prednaskach se jmenuje matNorm3)*)
Quiet@Remove["Global`*"];
n=150;
matSoustavy=Table[Prime[i+j-1],{i,1,n},{j,1,n+1}];
Export["matSoustavy.csv",matSoustavy];
matSoustavy=Import["matSoustavy.csv"];
n=Dimensions[matSoustavy][[1]];
A=N[Take[#,n]&/@matSoustavy];
b=N[Last/@matSoustavy];
sol=LinearSolve[A,b]//N;
matNorm3[mat_]:=Module[{},
Max[Abs[N[#]]&@Eigenvalues[mat]]];
cond[mat_,matNorm_]:=matNorm[Inverse@N@mat]*matNorm[N@mat];
{cond[A,matNorm3],Total[sol]}

(*V souboru "infoVykony.csv" jsou data ve tvaru {T1,T2,P} tak, ze plati ze pro Tϵ<T1,T2) odebirame vykon P v kW*)
(*v souboru "teplotyPulhodinove.csv" jsou prumerne pulhodinove teploty*)
(*urcete odebranou energii za dane obdobi*)
Quiet@Remove["Global`*"];
teploty=Flatten[Import["teplotyPulhodinove.csv"]];
fPom[T_]:=a+b*(1/(T+30));

solab=Solve[{fPom[13]==0,fPom[-20]==8.}][[1]];
P[T_]:=(fPom[T]/.solab)/;T≤13;
P[T_]:=0./;T>13;
Plot[P[T],{T,-20,20}];
myR:=0.1*Round[10#]&;
infoVykony=Table[{i,i+1,myR[P[i+0.5]]},{i,-30,50}];
Export["infoVykony.csv",infoVykony];
infoVykony=Import["infoVykony.csv"];
pom=infoVykony/.{T1_,T2_,P_}:>{P,T1≤T<T2};
ClearAll[vykon];
vykon[T_]=Piecewise[pom];
Total[vykon/@teploty]*1000*1800

(*v souboru "teplotyPulhodinove.csv" jsou prumerne pulhodinove teploty*)
(*Urcete poradi dni v tydnu(*prvni den znamena prvnich 48 hodnot pulhodinovych teplot*) podle vzajemne korelovanosti prumernych dennich teplot*)
(*tedy ve Wolframu Correlation[listPondelnichTeplot,listPondelnichTeplotOTydenPozdeji]*)
Quiet@Remove["Global`*"];
SetDirectory[NotebookDirectory[]];
teploty=Flatten[Import["teplotyPulhodinove.csv"]];
denniPrumery=Mean/@Partition[teploty,2*24];
korelace[den_Integer]:=Module[{list1,list2,nMax},
nMax=Length[denniPrumery];
list1=Table[denniPrumery[[i+den-1]],{i,1,nMax-13,7}];
list2=Table[denniPrumery[[i+den-1+7]],{i,1,nMax-13,7}];
Correlation[list1,list2]
];
{#,korelace[#]}&/@Range[7];
SortBy[%,Last]

(*seriova kombinace civky a rezistoru se v case t0 pripoji ke zdroji napeti*)
(*pro jaky cas pripnuti vznikne nemensi prechodny dej, tedy nejmensi integral z kvadratu rozdilu prubehu a ustaleneho probehu?"*)
Quiet@Remove["Global`*"];
L=0.1;
R=50;
f=50;
T=(1/f);
ω=2*Pi*f;
tmax=T;
uz[t_]:=230.*(Sqrt[2])*Sin[ω*t];
iHUS=Im[Exp[I*ω*t]*((230.*(Sqrt[2]))/(I*ω*L+R))];
pl1=Plot[iHUS,{t,0,tmax},AxesLabel->{"t[s]","i[A]"},PlotStyle->{Black,Thickness[0.006]},GridLines->Full];
tStart=t/.FindRoot[iHUS==0,{t,0.005}];
rce[t0_]:={i[t0]==0,R*i[t]+L*i'[t]==uz[t]};
ukaz[t0_]:=Module[{sol,pl},
sol=NDSolve[rce[t0],i,{t,t0,tmax}][[1]];
pl=Plot[i[t]/.sol,{t,t0,tmax},AxesLabel->{"t[s]","i[A]"},PlotStyle->{Red,Thickness[0.006],Dashed},GridLines->Full];
Show[pl1,pl]
];
ukaz[tStart]

(*seriova kombinace kondenzatoru a rezistoru se v case t0 pripoji ke zdroji napeti*)
(*pro jaky cas pripnuti vznikne nemensi prechodny dej, tedy nejmensi integral z kvadratu rozdilu prubehu a ustaleneho probehu?"*)
Quiet@Remove["Global`*"];
c=220*(10^(-6));
R=50;
f=50;
T=(1/f);
ω=2*Pi*f;
tmax=2T;
uz[t_]:=230.*(Sqrt[2])*Sin[ω*t];
iHUS=Im[Exp[I*ω*t]*((230.*(Sqrt[2]))/((1/(I*ω*c))+R))];
pl1=Plot[iHUS,{t,0,tmax},AxesLabel->{"t[s]","i[A]"},PlotStyle->{Black,Thickness[0.006]},GridLines->Full];
tStart=t/.FindRoot[iHUS==0,{t,0.005}];
rce[t0_]:={i[t0]==0,R*i[t]+uc[t]==uz[t],i[t]==c*uc'[t]};
ukaz[t0_]:=Module[{sol,pl},
sol=NDSolve[rce[t0],i,{t,t0,tmax},StartingStepSize->(10^(-6))][[1]];
pl=Plot[i[t]/.sol,{t,t0,tmax},AxesLabel->{"t[s]","i[A]"},PlotStyle->{Red,Thickness[0.006],Dashed},GridLines->Full];
Show[pl1,pl]];
ukaz[tStart]