function [Eopt, Esimple]= runSimulink_new_formulation2
filename=strcat(pwd,'\OPTdata\Dummy03-Jul-2013.mat');
aux=load(filename);
 sol=aux.sol;
 
finaltime=sol.t(end);
tf=sol.t(end);
pricesim.time= sol.t;
pricesim.signals.values=sol.price;
pricesim.signals.dimensions=[1];

pricemean=mean(sol.price);


qshower=sol.qout.*(sol.T-5-273)/(50-5);
qshowersim.time= sol.t;
qshowersim.signals.values=qshower;
qshowersim.signals.dimensions=[1];






Qsim.time= sol.t;
Qsim.signals.values=sol.Q;
Qsim.signals.dimensions=[1];





qinsim.time= sol.t;
qinsim.signals.values=sol.qin;
qqinsim.signals.dimensions=[1];


TminS=50;
TmaxS=90;
VminS=0.05;
VmaxS=0.2;
Vmax=0.2;
Vmin=0.05;
Tmin=TminS+273;
Tmax=TmaxS+273;
Qmax=5.5;
[x0,par]=getData;
cp=par.cp;
Tin=par.Tin;

rho=1000;
cp=4.19;
qshower=1*4.16;
twc=60*30;
hout=1.5*rho*cp*qshower*(50-5)/(60*1000);
Hlb=rho*Vmin*cp*(Tmin-Tin);
Hub=rho*Vmax*cp*(Tmax-Tin);


deltat=0.5*(Hub+Hlb)/(hout-Qmax)
%case II
% hout1=rho*cp*4.19*(50-5)/(60*1000);
% twc1=60*45;
% hout2=rho*cp*2*4.19*(50-5)/(60*1000);
% twc2=60*20;


%case I
hout1=rho*cp*1.5*4.16*(50-5)/(60*1000);
twc1=60*25;
% twc1=deltat;
hout2=rho*cp*1*1.5*4.16*(50-5)/(60*1000);
twc2=60*25;

tworst=20*60;
hworst=(Hub-Hlb)/(2*tworst)

Hmin=max(twc1*(hout1-Qmax)+Hlb,twc2*(hout2-Qmax)+Hlb)

% Hmin=tworst*(hworst-Qmax)+Hlb
Hcritical=(Hub-Hlb)*0.5+Hlb



Vlow=(Hcritical)/(cp*rho*(273+90-Tin));
Vhigh=(Hcritical)/(cp*rho*(273+50-Tin));
Vsp=(Vlow+Vhigh)/2;

% Loss_max=3600*(max(sol.price)-min(sol.price))*twc2*(hout2-Qmax)/3600


Hmax=rho*cp*Vmax*(Tmax-Tin)
Hmin=rho*cp*Vmin*(Tmin-Tin)
Thw=50+273;
longDist=4*3600;
Qdemandmax=Qmax;
qdemandmax=0.9*Qdemandmax/(rho*cp*(Thw-Tin))*(60*1000)

% qshower=[qdemandmax qdemandmax];
% peaktimes=[0 5  5.+2.2  16 16.85 26]*3600;
% qshowervec=[zeros(100,1) ; zeros(100,1)+qshower(1) ;   zeros(100,1) ;zeros(100,1)+qshower(2) ;zeros(100,1)   ];
% timevector=[linspace(peaktimes(1),peaktimes(2),100)'  
%             linspace(peaktimes(2),peaktimes(3),100)'
%             linspace(peaktimes(3),peaktimes(4),100)'
%             linspace(peaktimes(4),peaktimes(5),100)'
%             linspace(peaktimes(5),peaktimes(6),100)'];

% timeshower=sol.t(sol.qout>0);
% deltashower=timeshower(end)-timeshower(1);

%enery-out rate
%  hout=max(qshower/(60*1000))*(50-5)*rho*cp;
%  Hout=hout*deltashower/3600


% qshowersim.time= timevector;
% qshowersim.signals.values=qshowervec;
% qshowersim.signals.dimensions=[1];


  options = simset('SrcWorkspace','current');
  sim('closed_loop_new.slx',[],options)

% sim('closed_loop_new.slx')
%%
close all

time=timeSim/3600;
        VscaledOpt=(Vopt-VminS)/(VmaxS-VminS);
        VscaledSimple=(Vsimple-VminS)/(VmaxS-VminS);
        TscaledOpt=(Topt+273-Tin)/(TmaxS+273-Tin);       
        TscaledSimple=(Tsimple+273-Tin)/(TmaxS+273-Tin);    
        Tinscaled=(Tin-Tin)/(TmaxS-Tin);

        
        capopt=Vopt*cp.*(Topt+273-Tin);
        capsimple=Vsimple*cp.*(Tsimple+273-Tin);
        capgreedy=Vgreedy*cp.*(Tgreedy+273-Tin);
        
        capmax=cp*Vmax*(Tmax-Tin);
        capmin=cp*Vmin*(Tmin-Tin);
        
        
       
        capopt=(capopt-capmin)/(capmax-capmin);
        capsimple=(capsimple-capmin)/(capmax-capmin);
        capgreedy=(capgreedy-capmin)/(capmax-capmin);
        HminS=Hcritical/rho;
        HminS=(HminS-capmin)/(capmax-capmin);
%           capopt=VscaledOpt.*cp.*(TscaledOpt-Tinscaled);
%         capsimple=VscaledSimple.*cp.*(TscaledSimple-Tinscaled);
%         capmax=cp*1*(1)*ones(size(capopt));
        
        
    
       Loss_simple=(Jsimple(end)-Jopt(end))
        Loss_greedy=(Jgreedy(end)-Jopt(end))
        
        

        
% capopt=cp*Vopt.*(Topt+273-Tin);
% capsimple=cp*Vsimple.*(Tsimple+273-Tin);


xlimite=[0 28];
linewidth=1.5;
Colorgreedy=rgb('DarkGreen');
Colorgreedy=[0 0.5 0];

figure(1)
subplot(3,1,1)
plot(time,1000*Vgreedy,'Color',Colorgreedy,'LineWidth',linewidth)
hold on
plot(time,1000*Vopt,'r-',time,1000*Vsimple,'b','LineWidth',linewidth)
set(gca,'FontSize',14)
ylabel('$V$, $l$','interpreter','latex')
axis tight
% ylim([0.05 0.21])
ylim([45 210])
xlim(xlimite)
subplot(3,1,2)
plot(time,Tgreedy,'Color',Colorgreedy,'LineWidth',linewidth)
hold on
plot(time,Topt,'r-',time,Tsimple,'b','LineWidth',linewidth)
set(gca,'FontSize',14)
ylabel('$T$, $^{\circ}C$','interpreter','latex')
axis tight
ylim([TminS-3 TmaxS+3])
xlim(xlimite)
subplot(313)
plot(time,HminS*ones(size(time)),'--','Color',rgb('Gray'),'Linewidth',1)
hold on
plot(time,capgreedy,'Color',Colorgreedy,'LineWidth',linewidth)
hold on
plot(time,capopt,'r-',time,capsimple,'b','LineWidth',linewidth)
set(gca,'FontSize',14)
ylabel('$E, scaled$','interpreter','latex')
xlabel('Time, $h$','interpreter','latex')
 text(3,0.3,...
      '$E_{\textit{backoff}}$',...
      'FontSize',16,'interpreter','latex')
% text(3.25,HminS+0.12,...
%      '$\uparrow$',...
%      'FontSize',16,'interpreter','latex')
 
 
axis tight
ylim([0 1.05])
xlim(xlimite)
%   saveas(gcf, 'C:\Users\viniciud\Documents\Vinicius\PhD NTNU\my papers\NTNU\JPC\pics\states_vary_price', 'epsc') %Save figure



figure(2)
subplot(3,1,1)
plot(time,Qgreedy,'Color',Colorgreedy,'LineWidth',linewidth)
hold on
plot(time,Qopt,'r-',time,Qsimple,'b','LineWidth',linewidth)
set(gca,'FontSize',14)
ylabel('$Q$, $kW$','interpreter','latex')
axis tight
ylim([0  5.7])
xlim(xlimite)
subplot(3,1,2)
plot(time,qingreedy,'Color',Colorgreedy,'LineWidth',linewidth)
hold on
plot(time,qinopt,'r-',time,qsimple,'b','LineWidth',linewidth)
set(gca,'FontSize',14)
ylabel('$q_{in}$, $l/min$','interpreter','latex')
axis tight
xlim(xlimite)
subplot(3,1,3)
plot(time,Jgreedy,'Color',Colorgreedy,'LineWidth',linewidth)
hold on
plot(time,Jopt,'r-',time,Jsimple,'b','LineWidth',linewidth)
set(gca,'FontSize',14)
ylabel('$Cost$, $ \$ $','interpreter','latex')
xlabel('Time, $h$','interpreter','latex')
axis tight
xlim(xlimite)
%  ylim([0 70])
 

Eopt=Jopt(end);
Esimple=Jsimple(end);
 
%    saveas(gcf, 'C:\Users\viniciud\Documents\Vinicius\PhD NTNU\my papers\NTNU\JPC\pics\inputs_vary_price', 'epsc') %Save figure


figure(3)
% subplot(213)
% plot(time,objectivegreedy,'Color',Colorgreedy,'LineWidth',linewidth)
% hold on
% plot(time,objectiveOpt,'r-',time,objectiveSimple,'b','LineWidth',linewidth)
% set(gca,'FontSize',14)
% ylabel('$Objective$','interpreter','latex')
% xlabel('Time, $h$','interpreter','latex')
% ylim([0 60])
% xlim(xlimite)
subplot(212)
plot(time,prices*3600,'Color',rgb('Tan'),'LineWidth',linewidth)
set(gca,'FontSize',14)
ylabel('$Price$, $ \$/kWh $','interpreter','latex')
xlabel('Time, $h$','interpreter','latex')
% ylim([0.5 3.2])
xlim(xlimite)
 subplot(211)
 plot(time,qshower,'Color',rgb('Tan'),'LineWidth',linewidth)
 set(gca,'FontSize',14)
ylabel('$q_{hw}$, $l/min$','interpreter','latex')
% xlabel('Time, $h$','interpreter','latex')
ylim([0.0 9])
xlim(xlimite)
%  saveas(gcf, 'C:\Users\viniciud\Documents\Vinicius\PhD NTNU\my papers\NTNU\JPC\pics\dist_vary_price', 'epsc') %Save figure
 
tilefigs



EnergyGreedy=trapz(time,Qgreedy)
timePeak=time(time<10);
ESoc=trapz(timePeak,Qsimple(time<10))
EnergyOpt=trapz(timePeak,Qopt(time<10))
OnPeakReduction=(ESoc-EnergyGreedy)/EnergyGreedy
OnPeakReductionOpt=(EnergyOpt-EnergyGreedy)/EnergyGreedy

