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= 150;                      % maximální prùtok turbínou
Vmaxout = 150*60*60;                % Maximální hodinový odtok z pøehrady
P=72;                               % výkon celkový
l = 4000;                           % délka (wiki)
d=45;                              % ¹íøka (wiki)
ucinturb = 0.9;
ucinalt = 0.96;
ucin = ucinturb*ucinalt;
horizont = 24*13;                     % Èasový horizont 
Cv = round(rand(horizont,1)*10+30)';  % cena elektøiny 33-37 ojro
tan = tand(80);
V0 = tan*l*hp*hp + d*l*hp;
VMAX = tan*l*Hmax*Hmax + d*l*Hmax;
VMIN = tan*l*Hmin*Hmin + d*l*Hmin;
num = xlsread('/Users/Martin/Downloads/ceny.xlsx','A19705:D30648');
Cv = num(1:horizont,3);
Qin = num(1:horizont,4)';

% Matice cen elektøiny...Diagonální matice s prvky ceny v jednotlivých
% hodinách na trhu
for k=1:horizont
    C(k,k)= Cv(k);
end



%promìnné, H....vý¹ka, Q...prùtok, UC....unit
UC = binvar(1,horizont);
Q  = sdpvar(1,horizont);
H  = sdpvar(1,horizont);
S  = sdpvar(1,horizont);

% Objekt Constr = podmínky pro øe¹ení optimlizaèního problému

Constr = [];
% Podmínky pro urèení maximálního a minimálního prùtoku
Constr = [Constr, UC*Qmin <= Q <=UC*Qmax];
% Podmínky pro urèení maximální a minimální vý¹ky hladiny
Constr =[Constr, Hmin <= H <= Hmax];
% Podmínky pro urèení vlivu pøítokù na 
Constr =[Constr, Smin <= S <= Smax];
for k=1:horizont   
    
    if k>1
        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];
        
    else
%         Constr = [Constr, H(1,k)==hp];
        Constr = [Constr, l*(H(1,k)-hp)*(2*tan*hp+d) == Qin(1,k)*3600 - 3600*Q(1,k)-3600*S(1,k)];
    end
end

% Constr = [Constr, V(1,k)==tan*H(1,k)*H(1,k)*l + d*l*H(1,k)];
obj = -(Q*0.0098+0.0001)*C*H'*ucin;

% ops = sdpsettings('solver','bmibnb','bmibnb.maxiter',10,'bmibnb.uppersolver','fmincon');
ops = sdpsettings('solver','bmibnb','bmibnb.uppersolver','fmincon');
sol = optimize(Constr,obj,ops);
a= value(Q);
b= value(H);
c= value(V);

stairs(a)
hold on
stairs(b,'black')
hold on
stairs(Cv,'g')
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/MIQP_13den')
savefig('/Users/Martin/Desktop/Matlab_diplomka/Vysledky/MIQP_13den')
