%%
clear all

% Eingangswerte definieren einen Tag 
%Zulauf = [3 6 4 5 4 10 30 36 46 80 84 74 65 53 64 53 74 47 38 21 10 6 10 3]; % Input V_15
%Preis = [2 2 2 3 4 6 9 15 20 50 24 34 14 16 5 13 26 34 18 13 6 4 3 2]; %Preise für einen Tag

% Eingangsdaten eine Woche
load('Eingang')
Zulauf=round(zl);
Preis= Preise;

nPeriods=length(Zulauf);            % Anzahl der Zeitschritte
nPump=2;                % Anzahl der Pumpen
Volumenstrom=[50;70];   % Volumenstrom der zwei Pumpem
c=0.1;                  % Verhältniss zwischen Volumenstrom und Leistungsaufnahme
Last=Volumenstrom*c;    % Verbrauch der Pumpen

%% Zielfunktion/ Kostenfunktion f =[73x1]

f_pump1= Preis'*Last(1);      % Preis Pumpe1 Durchsatz und cp Werte Durchfluss-> Leistung
f_pump2= Preis'*c*Last(2);
f_fs = zeros(nPeriods+1,1);             % an den Füllstand sind keine Kosten
                                        % geknüpft Startwert für Füllstand muss
                                        % jedoch integriert werden

f=[f_pump1;f_pump2;f_fs];               % Zusammensetzen des Vektors für die Kostenfunktion

%% Bounds

lb_p=zeros(nPeriods,nPump);     % Pumpen ist aus
ub_p=ones(size(lb_p));          % Pumpen sind an
lb_fs=ones(nPeriods+1,1)*10;    % Füllstand muss zwischen 10 und 400 liegen
ub_fs=ones(nPeriods+1,1)*400;
ub_fs(end)=200;
lb=[lb_p(:);lb_fs(:)];          % einen Vektor erzeugen für gesamt X Vektor
ub=[ub_p(:);ub_fs(:)];

%% A und b definieren 

Startwert= 150; % Füllstand bei dem gestartet wird.
b=[Startwert;Zulauf'];

% Matrix A Formulieren mit A=[A1 A2 A3] 

% A1= Volumenstrom welcher abgeführt werden kann wenn x(t)= 1 ist
% A2= Volumenstrom der abgeführt werden kann wenn x(t)= 1 ist
% A3= Koeffizieten der Füllstände 1 oder -1 

A1=[zeros(1,nPeriods);eye(nPeriods)*Volumenstrom(1)];
A2=[zeros(1,nPeriods);eye(nPeriods)*Volumenstrom(2)];
A3=eye(nPeriods+1);

for i=2:nPeriods+1
    A3(i,i-1)=-1;
end
A=[A1 A2 A3];
A=sparse(A); % Nullen aus der matrix entfernen
%% Binäre Variablen Festlegen 

intcon= 1:nPeriods*2;

%% Gleichungssystem lösen
options=optimoptions(@intlinprog,'TolGapRel',1e-2,'MaxNumFeasPoints',4);

[x,Kosten]= intlinprog(f,intcon,[],[],A,b,lb,ub,options)
%[x,Kosten]= intlinprog(f,intcon,[],[],A,b,lb,ub)


%% Werte darstellen
% 
% Pumpen(:,1)=x(1:nPeriods);
% Pumpen(:,2)=x(nPeriods+1:nPeriods*2);
% fs=x(nPeriods*2+1:end);
% 
% %Pumpemaktivität 
% figure(1)
% subplot(2,1,1)
% bar(1:nPeriods,Pumpen)
% title('Aktivität der Pumpen')
% ylim([0 1.2]);
% grid on
% 
% l=legend({'P1 klein','P2 groß'},'FontSize',8,'FontWeight','bold','Location', 'NorthWest');
% 
% subplot(2,1,2)
% H1=area(1:nPeriods,Pumpen(:,1)*Volumenstrom(1));
% grid on
%  set(H1,'FaceColor','r');
% set(get(H1,'Children'),'FaceAlpha',0.8);
% hold on
% H2=area(1:nPeriods,Pumpen(:,2)*Volumenstrom(2));
% set(get(H2,'Children'),'FaceAlpha',0.8)
% 
% 
% 
% ph(:,1)=Pumpen(:,1)*Volumenstrom(1);
% ph(:,2)=Pumpen(:,2)*Volumenstrom(2);
% 
% figure(3)
% area(1:nPeriods,ph)
% 
% % Speicherstand vgl Preise
% figure(2)
% 
%  [AX,H1,H2]=plotyy(1:nPeriods,fs(2:end),1:nPeriods,Preis);
%  grid on
% set(get(AX(1),'Ylabel'),'String','Speicherstand')
% set(get(AX(1),'Ylabel'),'FontSize',14)
% set(get(AX(2),'Ylabel'),'String','Preise')
% set(get(AX(2),'Ylabel'),'FontSize',14)
% 
% xlabel('Stunden/Woche')
% title('Vergleich Speicherstand/ Pumpenaktivität ','FontSize',14)
% 
%  set(H1,'LineWidth',1,'Color',[1 0 0]);
% % set(get(H1,'Children'),'FaceAlpha',0.2);
% set(H2,'LineWidth',1,'Color',[0 0 1]);
% set(AX(1),'FontSize',14,'YColor',[1 0 0]);
% set(AX(2),'FontSize',14,'YColor',[0 0 1]); 
% set(AX(1),'YTick',[0,50,100,200,300,400,500]);
% 
%  l=legend({'Speicherstand','Preisverlauf'},'FontSize',8,'FontWeight','bold','Location', 'NorthWest');
% 
%