%Main file to simulate room temperature
%October 2018
%This simulates the temperataure of a room where:
%MVs:
%QH1= heated water [kJ/s]
%QH2= electric heating [kJ/s]
%QAC= air conditioning [kJ/s]


%Disturbance:
%Tin= outdoor temperature that goes through ventilation

%States:
%T= indoor temperature; x1
%Tfl= floor temperature; x2

%The model is in the file called "Room". 
%The Simulink file is called RoomT
%Open loop Simulink file is called RoomT_OL

clc
clear
close all

%% Nominal values
QH10    =0;         %at steady-state, no need of heating water
QH20    =0;         %at steady-state, no need of electric heating
QAC0    =0;         %at steady-state, no need of air conditioning
m_in0   =0.4;       %mass flow for ventilation; kg/s

Qhfl0=0;            %at steady-state, no heat transfer through the floor

%Nominal values for states
T0      =21;          %room temperature, °C
Tfl0    =21;        %floor temperature, °C
x0      =[T0; Tfl0]; 

%Disturbance: ambient temperature
Tin0    = 21;        %°C / going from -10 to 30°C

%% Max and min val
HWmin   = 0;
HWmax   = 3;

EHmin   = 0;
EHmax   = 4;

ACmin   = 0;
ACmax   = 4.5;

mvmax=0.4;
mvmin=0.1;
mv_sp=0.4;


%% Parameters 
% Parameters for model
m       =200;          %mass of air in the room [kg]
Cp      =1;           %Cp of air [kJ/(kg K)]
Ufl     =10/1000;    %heat transfer coefficient for floor [kW/m^2K]=10 W/m^2K
Afl     =50;         %area of floor [m^2]
mfl     =600;        %mass of floor [kg]
Cpfl    =1;         %Cp of floor - concrete [kJ/(kg k] 

%Parameters for objective function
% actual electricity prices here: https://www.nordpoolgroup.com/
price   = 500;

p_el    = 1; %price/1000;  %    % 290 kr/MWh
p_hw    = 2/3; %0.5*price/(1000*0.9);  %    1.5; %0.05; 
p_ac    = 1/3; %price/(1000*3);  %Assume energy efficiency of 3
kT      = 0.2; %0.05;
kv      = 0;
%Parameter vectors for simulations
p       = [m,Cp, Ufl, Afl, mfl, Cpfl]';
par_el  = [p_el, m_in0, Cp, kT]';
par_hw  = [p_hw, m_in0, Cp, kT]';
par_ac  = [p_ac, m_in0, Cp, kT]';
%par_mv  = [kv, kT, Cp, EHmax, ACmax,m_in0]';
par_J   = [p_el, p_hw, p_ac, m_in0, kv, kT]';



%% Desired temperature
Tdes0 =  T0;
%% Steps in disturbance Tin
Tin_1   = 10;
Tin_2   = -5;
Tin_3   = -8;
Tin_4   = -13;
Tin_5   = 0;

tTin_1  = 3600/2;
tTin_2  = 3.5*3600;
tTin_3  = 7*3600;
tTin_4  = 10*3600;
tTin_5  = 20000;

%% Room Tuning
RoomTuning

%% SRC
StdSRC_3MVs
%% Simulate
tend=15*3600; %max([tQH1_2,tQAC_2,tQH2_2,tm_in_2,tTin_2,tm_in_1,tm_in_2])+50000;  %simulation time
sim('RoomT_SRC_gen')

%% Plot
%t is in seconds
tplot=t/3600;
timescale='h';

ACmaxplot = ones(length(tplot),1)*ACmax;
HWmaxplot = ones(length(tplot),1)*HWmax;
EHmaxplot = ones(length(tplot),1)*EHmax;
mvmaxplot = ones(length(tplot),1)*mvmax;

Qminplot = zeros(length(tplot),1);
mvminplot = ones(length(tplot),1)*mvmin;

set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',24,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',24,...
'DefaultLineLineWidth',2.5,...
'DefaultLineMarkerSize',12)

figure
subplot(3,1,1)

plot(tplot,Tdes_sp,'--','Linewidth', 1,'Color', [0.7500         0    0.7500] )
hold on
plot(tplot,T,'Color', [ 0.4940    0.1840    0.5560])
ylabel('T [°C]')
ylim([15 24])
title ('$y = T$','Interpreter','Latex');
legend('$T^{des}$','$T$','Interpreter','Latex', 'Orientation','horizontal');
hold off
legend boxoff 

subplot(3,1,2)
plot(tplot,Tdes_sp,'--','Linewidth', 1,'Color', [ 0.7500         0    0.7500] )
hold on
plot(tplot,Tin, 'Color', [ 0.4660    0.6740    0.1880])
ylabel('T, [°C]')
ylim([0 35])
title ('$d = T^{amb}$','Interpreter','Latex');
legend('$T^{des}$','$T_{in}$','Interpreter','Latex','Orientation','horizontal');
legend boxoff 

subplot(3,1,3)
plot(tplot,T, 'Color', [ 0.4940    0.1840    0.5560])
hold on
plot(tplot,Tfl, '-.','Color', [   0    0.7500    0.7500])
plot(tplot,Tdes_sp,'--','Linewidth', 1,'Color', [ 0.7500         0    0.7500] )
title ('$x = [T, T_{fl}]$','Interpreter','Latex');
ylabel('T, [°C]')
legend('$T$', '$T_{fl}$','$T^{des}$','Interpreter','Latex','Orientation','horizontal');
hold off
legend boxoff 
ylim([15 30])
xlabel(strcat('time [',timescale,']'))

%error
 figure
plot(tplot,e)
hold on
plot(tplot,zeros(length(t),1))
%%
%sugggested MVs
figure
plot(tplot,uAC1)
hold on
plot(tplot,uHW1)
plot(tplot,uEH1)
legend('AC','HW', 'EH')
%%
%MVs
figure
subplot (3,1,1)
p1=plot(tplot,uAC);
hold on
plot(tplot,ACmaxplot,'--', 'Color',[0.5 0.5 0.5],'Linewidth', 1)
plot(tplot,Qminplot,'--', 'Color',[0.5 0.5 0.5],'Linewidth', 1)
ylabel('$\dot{Q}_{ac}$ [kW]','Interpreter','Latex');
title('$\dot{Q}_{ac}$','Interpreter','Latex');
ylim([-0.2 ACmax+0.2])
ytickformat('%.1f')
%legend('$\dot{Q}_{ac}$','$\dot{Q}_{ac}^{max}$','Interpreter','Latex');
%legend([p1],'$\dot{Q}_{ac}$','Interpreter','Latex');
hold off
%legend boxoff 

subplot (3,1,2)
plot(tplot,uHW,'Color', [ 0.9290    0.6940    0.1250] )
hold on
plot(tplot,HWmaxplot,'--', 'Color',[0.5 0.5 0.5],'Linewidth', 1)
plot(tplot,Qminplot,'--', 'Color',[0.5 0.5 0.5],'Linewidth', 1)
ylabel('$\dot{Q}_{hw}$ [kW]','Interpreter','Latex');
title('$\dot{Q}_{hw}$','Interpreter','Latex');
ylim([-0.2 HWmax+0.2])
ytickformat('%.1f')
%title('Heating water')
%legend('$\dot{Q}_{hw}$','$\dot{Q}_{hw}^{max}$','Interpreter','Latex');
hold off
%legend boxoff 

subplot (3,1,3)
plot(tplot,uEH , 'Color', [1.0000         0         0])
hold on
plot(tplot,EHmaxplot,'--', 'Color',[0.5 0.5 0.5],'Linewidth', 1)
plot(tplot,Qminplot,'--', 'Color',[0.5 0.5 0.5],'Linewidth', 1)
ylabel('$\dot{Q}_{eh}$ [kW]','Interpreter','Latex');
title('$\dot{Q}_{eh}$','Interpreter','Latex');
ylim([-0.2 EHmax+0.2])
ytickformat('%.1f')
%title('Electrical heating')
%legend('$\dot{Q}_{el}$','$\dot{Q}_{el}$','Interpreter','Latex');
hold off
%legend boxoff 
xlabel(strcat('time [',timescale,']'))
%%


figure
plot(tplot, Jactual)
title('Objective function')

% %% figure
% subplot(2,1,1)
% plot(tplot,T)
% hold on
% plot(tplot,Tfl)
% ylabel('Temperature, [°C]')
% legend('$T$', '$T_{fl}$','Interpreter','Latex');
% hold off
% legend boxoff 
% 
% subplot(2,1,2)
% plot(tplot,Tin)
% ylabel('Temperature, [°C]')
% xlabel('time, [min]')
% legend('$T_{in}$','Interpreter','Latex'); 
% legend boxoff 
% %%
% figure
% plot(tplot,ACsp,'--', 'Linewidth', 1)
% hold on
% plot(tplot,EHsp,'--', 'Linewidth', 1)
% plot(tplot,HWsp,'--', 'Linewidth', 1)
% plot(tplot,T)
% ylabel('T [°C]')
% %ylim([20.5 22])
% legend('$Q_{ac}^{sp}$', '$Q_{el}^{sp}$','$Q_{hw}^{sp}$', 'T','Interpreter','Latex');
% hold off
% legend boxoff 
% 
% %%
% figure
% subplot(2,1,1)
% plot(tplot,ACsp)
% hold on
% plot(tplot,EHsp)
% plot(tplot,HWsp)
% ylabel('T [°C]')
% ylim([20.5 22])
% legend('$Q_{ac}^{sp}$', '$Q_{el}^{sp}$','$Q_{hw}^{sp}$','Interpreter','Latex');
% hold off
% legend boxoff 
% 
% subplot(2,1,2)
% plot(tplot,m_v_sp)
% ylabel('mass flow, [kg/s]')
% xlabel('time, [min]')
% legend('$m_{v}^{sp}$','Interpreter','Latex');
% legend boxoff 
% 
% figure
% plot(tplot,uAC)
% hold on
% plot(tplot,uEH)
% plot(tplot,uHW)
% ylabel('MV [kW]')
% %ylim([20.5 22])
% legend('$Q_{ac}$', '$Q_{el}$','$Q_{hw}$','Interpreter','Latex');
% hold off
% legend boxoff 
% % 
% % figure
% % subplot(2,2,1)
% % plot(tplot,T)
% % hold on
% % plot(tplot,Tfl)
% % ylabel('Temperature, [°C]')
% % legend('$T_{room}$', '$T_{floor}$','Interpreter','Latex');
% % hold off
% % 
% % subplot(2,2,3)
% % plot(tplot,QAC,':')
% % hold on
% % plot(tplot,QH2)
% % plot(tplot,QH1,'--')
% % ylabel('Heat, [kW]')
% % legend('$Q_{AC}$', '$Q_{EL}$','$Q_{HW}$','Interpreter','Latex');
% % hold off
% % xlabel('time, [min]')
% 
% subplot(2,2,4)
% plot(tplot,m_in)
% ylabel('mass flow, [kg/s]')
% 
% legend('$m_{in}$','Interpreter','Latex');
% 
% subplot(2,2,2)
% plot(tplot,Tin)
% ylabel('Temperature, [°C]')
% xlabel('time, [min]')
% legend('$T^{amb}$','Interpreter','Latex'); 

%% Save results
% filename = 'RoomSRC.mat';
% save(filename)
% x = load('RoomSRC.mat');
% 
% names = fieldnames(x);
% for iname = 1:length(names)
%   x.(['RoomSRC_', names{iname}]) = x.(names{iname});
%   x = rmfield(x, names{iname});
% end
% 
% save('RoomSRC_.mat', '-struct', 'x');
