clc; clear all; Qin = 35.6; % Průměrný vtok m3/s Vinh = Qin*60*60; % Objem přítoku za hodinu Qmax = 150; % maximální průtok m3/s Qmin = 50; V=18.15E6; % objem nádrže v m3 Vmax = 39E6; % max objem nádrže dno = 45; % Šířka dna hp = 24.1; % průměrný spád na turbíny (wiki) Hmax = 25.8; Hmin = 23.4; Smin = 0; Smax = 1000; ro = 1000; g = 9.8; hvar = 1.7; % kolísání vodní hladiny Qmaxturb= 60; % maximální průtok turbínou Vmaxout = 150*60*60; % Maximální hodinový odtok z přehrady P=72; % výkon celkový ucinturb = 0.9; ucinalt = 0.96; ucin = ucinturb*ucinalt; l = 4000; % délka (wiki) d=45; % šířka (wiki) horizont = 24*365; % Časový horizont tan = tand(80); num = xlsread('/Users/Martin/Downloads/ceny.xlsx','A19705:D30648'); C = num(1:horizont,3)'; Qin = num(1:horizont,4)'; % generování srážek % proměnné, v tomto případě zanedbám proménnou h (jen LP) UC = binvar(1,horizont); Q = sdpvar(1,horizont); H = sdpvar(1,horizont); S = sdpvar(1,horizont); % Podmínky Constr = []; for k = 1:horizont Constr = [Constr, UC(1,k)*Qmin <= Q(1,k) <= UC(1,k)*Qmax]; Constr = [Constr, Hmin <= H(1,k) <= Hmax]; Constr = [Constr, Smin <= S(1,k) <= Smax]; end for k=1:horizont if k==1 Constr = [Constr, l*(H(1,k)-hp)*(2*tan*hp+d) == Qin(1,k)*3600 - 3600*Q(1,k)-3600*S(1,k)]; else Constr = [Constr, l*(H(1,k)-H(1,k-1))*(2*tan*hp+d) == Qin(1,k)*3600 - 3600*Q(1,k)-3600*S(1,k)]; Constr = [Constr, Q(1,k) - Q(1,k-1)<=80]; Constr = [Constr, Q(1,k-1) - Q(1,k)<=80]; end end Objective = 0; for k = 1:horizont Objective = Objective + ucin*C(k)*hp*(0.0098*Q(1,k)+ 0.001); end options = sdpsettings('verbose',1,'solver','gurobi'); sol = optimize(Constr,-1*Objective,options); value(H); value(Q); axis([0 2*pi -1.5 1.5]) a= value(Q); b=value(H); c= value(C); d= value(Qin); stairs(a) hold on stairs(b,'black') hold on stairs(c,'g') % hold on % stairs(d,'red') hold off axis([0 horizont -10 200]) xlabel('Časový horizont [h]') ylabel('Průtok [m3/s], Cena elektřiny [EUR], Výška hladiny [m]') legend('Průtok turbínou','Výška hladiny', 'Cena elektřiny') save('/Users/Martin/Desktop/Matlab_diplomka/Vysledky/MILP_rok') savefig('/Users/Martin/Desktop/Matlab_diplomka/Vysledky/MILP_rok') % stairs(Qin) % axis([0 horizont -10 150]) % xlabel('Časový horizont [h]') % ylabel('Přítoky[m3/s]')