
function [ J ] = objective_generall( decisionvector,start,days,price,time_excel)

%Finds the objective function
%Takes in the Tbufferi and Pthreshold, runs the simulink model and returns the cost

Tbuffervector = decisionvector(1);
Pthreshold_vector = decisionvector(2:end);


%Parameters
par.cp.water = 4.19; %[kJ/Kg,K] SI Chemical Data
par.ro.water = 1; %[g/cm^3] SI Chemical Data

%Saturation limits
Qmin = 0.0;     %[kW]
Qmax = 5.5;     %[kW]
qin_min = 0;    %[m^3/s]
qin_max = 10;   %[m^3/s]

Qlb = Qmin;
Qub = Qmax;
qinlb = qin_min;
qinub = qin_max;

load day.mat

            qhw=[];
            p=[];

for i = start:(start+days-1);
    
    tsim = 3600*24;
    
    %constant demand
        
    qhw = day.flow;
    Thws = zeros(size(day.time))'+273+50;
    Tin = zeros(size(day.time))'+273+5;
    
    
%     %constant price
%     day(i) = getProfile;   
%     qhw(:,i) = day(i).flow;
%     Thws(:,i) = zeros(size(day(i).time))'+273+50;
%     Tin(:,i) = zeros(size(day(i).time))'+273+5;
    
    
    
    % varying price
    % PRICE: converting excel values into the dim of the distrubance matrix
    
    price_day(:,i) = price(:,i);
    price_days(:,i) = interp1(time_excel,price_day(:,i),day.time,'linear','extrap');
    
    %constant price
    %p(:,i) = price_days(:,i)';
    
    %constant demand
    p=[p price_days(:,i)];
    
    %% Distubance matrix
%     
%     %constant price
%     dsim.time = [day(i).time]'*3600;
%     dsim.signals.values = [qhw(:,i) Thws(:,i) Tin(:,i) p(:,i)]; % disturbances
%     dsim.dimentions = [4 1];
    

    %constant demand
    dsim.time = [day.time]'*3600;
    dsim.signals.values = [qhw Thws Tin p]; % disturbances
    dsim.dimentions = [4 1];
     

    Tbuffertime = [0 5 9 14 17 20 21 24]'; %Time in h,x7, initial
    %Tbuffertime = [0 4 6 8 11 16 20 24]'; %Time in h,x71
    %Tbuffertime = [0 2 5 10 14 19 22 24]'; %Time in h, x72
    %Tbuffertime = [0 6 9 12 16 20 23 24]'; %Time in h, x73
    
    
    %case 5
    Tbuffer_vect = (zeros(size(Tbuffertime))+Tbuffervector')';
    Pthreshold_vect = (zeros(size(Tbuffertime))+Pthreshold_vector')';
    
    
    Pthreshold.time = [Tbuffertime]*3600;
    Pthreshold.signals.values = [Pthreshold_vect]';
    Pthreshold.dimentions = [1 1];
    
    
    Tbuffer.time = [Tbuffertime]*3600;
    Tbuffer.signals.values = [Tbuffer_vect]';
    Tbuffer.dimentions = [1 1];


    options = simset('SrcWorkspace','current');
    sim('simulinkmodel_Tbuffer',[],options)
    J=simout.signals.values(end,end);
    
    penalty = 100;
  
    numberofint = length(Tbuffervector);
    
     for m=1:numberofint;
      
    Jp = max(Tbuffervector(m)-90-273,0)*penalty+max(25+273-Tbuffervector(m),0)*penalty;
    
    J = J + Jp;
      
    end
   
    
    J;
    
    plotresPthres.qhw_case=simout1.signals.values;
    plotresPthres.JJ = simout.signals.values(:,3); 
    plotresPthres.T=simout.signals.values(:,2);
    plotresPthres.t = simout.time/3600;
    plotresPthres.Pthres = simout2.signals.values;
    plotresPthres.Tset = simout3.signals.values;
    
    save plotresPthres
   
end
end