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

%Finds the objective function
%Takes the proposed inital values and runs the simulink model 
%and returns the cost

%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 %Demand data

            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];
     
    %Tbuffer values
   
    Tbuffertime = [0 5 9 14 17 20 21 24]'; %Time in h,x7, initial
    
    %%case with different Time_int
    %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
   
    
   
    Tbuffer_vect = (zeros(size(Tbuffertime))+Tbufferi')';
    Pthreshold_vect = (zeros(size(Tbuffertime))+Pthreshold')';
    
    
    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;
    
    %Implemented penalty constraint
    
    numberofint = length(Tbufferi);
    
     for m=1:numberofint;
      
    Jp = max(Tbufferi(m)-90-273,0)*penalty+max(25+273-Tbufferi(m),0)*penalty;
    
    J = J + Jp;
      
     end
     
    T = simout.signals.values(end,2);
     
    plotres.qhw_case=simout1.signals.values;
    plotres.JJ = simout.signals.values(:,3); 
    plotres.T=simout.signals.values(:,2);
    plotres.t = simout.time/3600;
    plotres.Tset = simout2.signals.values;
    plotres.p = simout3.signals.values;
    
    %Case 4
    %Delta E penalty implemented for lacing final constraint.
    
%      deltaT = (90+273-T);
%      deltaE = par.cp.water*par.ro.water*deltaT*200*2.778*10^(-4)*0.127;
%     
%      J = J + deltaE;
%      
%      plotresJan.dEJJ =simout.signals.values(:,3) + deltaE;
   
save plotres
     
     J;
       
end
end