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]')

