%Optimal Temperature Control of rooms for Minimum Energy Cost

%Master thesis 2013

close all 
clear all
clc

% Function for finding the solution of optimization problem when using
% PI controllers
% Storage of heat in the floor
% Written by Siri Hofstad Trapnes, 2013

%% Parameters

%UA = J/hK
%mCp = J/K

par.UAfr = 0.350; % [kJ/sK = kW/K]
par.UAro = 0.007; % [kJ/sK = kW/K]
par.mCpf = 4000; % [kJ/K]
par.mCpr = 70; % [kJ/K]
par.mr = 70; % [kg]

min = 0.06; % [kg] 
mout = min;

%% Energy price

%Varying price
ts = 600; % 10min
tend = 20*60*60; % 24h
t = [0:ts:tend]'; % [s]
t_p = t/(60*60);

t1 = length(t(1):ts:t(37)); % start - 6h %
t2 = length(t(38):ts:t(49)); % 6 - 8 %
t3 = length(t(50):ts:t(end)); % 8 - end %

P1 = 1/(60*60);
P2 = 0.6/(60*60);

P = [ones(t1,1)*P2; ones(t2,1)*P1; ones(t3,1)*P2]; % 15h
P_plot = P;

var.time = t;  
var.signals.values = P;

%% Inputs and initial conditions

%Inputs
%qhr = 2; % No heat on: 0 %max: 2 [kW] 
%qhf = 2.5; % No heat on: 0 %max: 2.5 [kW] 
To = 273+0; % [K]

%Varying inputs
t1 = length(t(1):ts:t(25)); % start - 4h
t2 = length(t(26):ts:t(37)); % 4h - 6h
t3 = length(t(38):ts:t(end)); %15h

qhf_min = 0;
qhf_max = 2.5;
qhr_min = 0;
qhr_max = 2; 

qhf = [ones(t1,1)*qhf_min; ones(t2,1)*qhf_max; ones(t3,1)*qhf_min];
qhr = [ones(t1,1)*qhr_min; ones(t2,1)*qhr_max; ones(t3,1)*qhr_min];

varqhf.time = t;
varqhf.signals.values = qhf;

varqhr.time = t;
varqhr.signals.values = qhr;

%Initial conditions
Tf = 293; % [K] 
Tr = 292; % [K]
J = 0; % [NOK]
init = [Tf Tr J];

%Set point
Tf_set = 292; % [K]
Tr_set = 292; % [K]

%normal
Tfmax = 298; % [K]
Tfmin = 293; % [K] 293
Trmax = 298; % [K]
Trmin = 292; % [K] 293

%% Optimizing, storing of heat in the floor

%Optimal solution
%High price interval between 6 and 8

sim('simulink_master_PI_2')

%Outputs
Tf = simout(:,1)-273; % [C]
Tr = simout(:,2)-273; % [C]
J = simout(:,3); 
Qhf = qhf_out;
Qhr = qhr_out;
t = time/(60*60); % [h]

%Energy consumption
Jend = J(end)

%Solution from switching rule
ts = 600; % 10min
tend = 20*60*60; % 20h
t_switch = [0:ts:tend]'; % [s]
t_p_switch = t_switch/(60*60);

t1 = length(t_switch(1):ts:t_switch(37)); % start - 6h
t2 = length(t_switch(38):ts:t_switch(49)); % 6 - 8
t3 = length(t_switch(50):ts:t_switch(end)); % 8 - end

P = [ones(t1,1)*P2; ones(t2,1)*P1; ones(t3,1)*P2];
P_switch = P;
var.time = t_switch;  
var.signals.values = P;

sim('simulink_master_PI_2_switch')

%Outputs
Tf_switch = simout(:,1)-273; % [C]
Tr_switch = simout(:,2)-273; % [C]
J_switch = simout(:,3); 
Qhf_switch = qhf_out;
Qhr_switch = qhr_out;
t_switch = time/(60*60); % [h]
%Energy consumption
Jend_switch = J_switch(end)

% Broad peak
% High energy price from 5 h - 12 h

ts = 600; % 10min
tend = 20*60*60; % 20h
t_bp = [0:ts:tend]'; % [s]
t_p_broadpeak = t_bp/(60*60);

t1 = length(t_bp(1):ts:t_bp(31)); % start - 5h
t2 = length(t_bp(32):ts:t_bp(73)); % 5 - 12
t3 = length(t_bp(74):ts:t_bp(end)); % 12 - end

P = [ones(t1,1)*P2; ones(t2,1)*P1; ones(t3,1)*P2];
P_broadpeak = P;
var.time = t_bp;  
var.signals.values = P;

sim('simulink_master_PI_broadpeak')

%Outputs
Tf_broadpeak = simout(:,1)-273; % [C]
Tr_broadpeak = simout(:,2)-273; % [C]
J_broadpeak = simout(:,3); 
Qhf_broadpeak = qhf_out;
Qhr_broadpeak = qhr_out;
t_broadpeak = time/(60*60); % [h]

%Energy consumption
Jend_broadpeak = J_broadpeak(end)

%Solution from broad peak using switching rule

ts = 600; % 10min
tend = 20*60*60; % 24h
t_bp_s = [0:ts:tend]'; % [s]
t_p_broadpeak_s = t_bp_s/(60*60);

t1 = length(t_bp_s(1):ts:t_bp_s(31)); % start - 5h
t2 = length(t_bp_s(32):ts:t_bp_s(73)); % 5 - 12
t3 = length(t_bp_s(74):ts:t_bp_s(end)); % 12 - end

P = [ones(t1,1)*P2; ones(t2,1)*P1; ones(t3,1)*P2];
P_broadpeak_s = P;
var.time = t_bp_s;  
var.signals.values = P;

sim('simulink_master_PI_broadpeak_switch')

%Outputs
Tf_broadpeak_s = simout(:,1)-273; % [C]
Tr_broadpeak_s = simout(:,2)-273; % [C]
J_broadpeak_s = simout(:,3); 
Qhf_broadpeak_s = qhf_out;
Qhr_broadpeak_s = qhr_out;
t_broadpeak_s = time/(60*60); % [h]

%Energy consumption
Jend_broadpeak_s = J_broadpeak_s(end)

% %Broad peak 2. No time to heat up the floor
% %High energy price from 3 h - 10 h
% 
ts = 600; % 10min
tend = 20*60*60; % 20h
t_bp2 = [0:ts:tend]'; % [s]
t_p_broadpeak2 = t_bp2/(60*60);

t1 = length(t_bp2(1):ts:t_bp2(19)); % start - 3h
t2 = length(t_bp2(20):ts:t_bp2(61)); % 3 - 10
t3 = length(t_bp2(62):ts:t_bp2(end)); % 10 - end

P = [ones(t1,1)*P2; ones(t2,1)*P1; ones(t3,1)*P2];
P_broadpeak2 = P;
var.signals.values = P;

sim('simulink_master_PI_broadpeak2')

%Outputs
Tf_broadpeak2 = simout(:,1)-273; % [C]
Tr_broadpeak2 = simout(:,2)-273; % [C]
J_broadpeak2 = simout(:,3); 
Qhf_broadpeak2 = qhf_out;
Qhr_broadpeak2 = qhr_out;
t_broadpeak2 = time/(60*60); % [h]

%Energy consumption
Jend_broadpeak2 = J_broadpeak2(end)

% %Solution from broad peak 2 using switching rule
% 
% ts = 600; % 10min
% tend = 20*60*60; % 20h
% t_bp2_s = [0:ts:tend]'; % [s]
% t_p_broadpeak2_s = t_bp2_s/(60*60);
% 
% t1 = length(t_bp2_s(1):ts:t_bp2_s(19)); % start - 3h
% t2 = length(t_bp2_s(20):ts:t_bp2_s(61)); % 3 - 10
% t3 = length(t_bp2_s(62):ts:t_bp2_s(end)); % 10 - end
% 
% P = [ones(t1,1)*P2; ones(t2,1)*P1; ones(t3,1)*P2];
% P_broadpeak2_s = P;
% var.time = t_bp2_s;  
% var.signals.values = P;
% 
% sim('simulink_master_PI_broadpeak2_switch')
% 
% %Outputs
% Tf_broadpeak2_s = simout(:,1)-273; % [C]
% Tr_broadpeak2_s = simout(:,2)-273; % [C]
% J_broadpeak2_s = simout(:,3); 
% Qhf_broadpeak2_s = qhf_out;
% Qhr_broadpeak2_s = qhr_out;
% t_broadpeak2_s = time/(60*60); % [h]
% 
% %Energy consumption
% Jend_broadpeak2_s = J_broadpeak2_s(end)

%% Plots

%Plot of optimal solution and solution from switching rule
figure(1)
subplot(5,1,1)
stairs(t_p,P_plot*60*60,'b')
ylabel('Energy price [kr/kWh]')
xlabel('Time [h]')
title('Energy price, P') 

subplot(5,1,2)
hold on
plot(t,Tf,'b')
plot(t_switch,Tf_switch,'r')
plot(t,ones(length(t),1)*Tfmax-273,'k --')
plot(t,ones(length(t),1)*Tfmin-273,'k --')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tf')
legend('Optimal solution','Switching rule','T_{max}','T_{min}')
hold off

subplot(5,1,3)
hold on
plot(t,Qhf,'b')
plot(t_switch,Qhf_switch,'r')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhf') 
hold off

subplot(5,1,4)
hold on
plot(t,Tr,'b')
plot(t_switch,Tr_switch,'r')
plot(t,ones(length(t),1)*Trmax-273,'k --')
plot(t,ones(length(t),1)*Trmin-273,'k --')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tr')  
hold off

subplot(5,1,5)
hold on
plot(t,Qhr,'b')
plot(t_switch,Qhr_switch,'r')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhr')
hold off

%Saving in tikz
%matlab2tikz('minimumPointsDistance',{0.01},'optimal_solution_TfTruP.tikz')

figure(2)
hold on
plot(t,J,'b')
plot(t_switch,J_switch,'r')
ylabel('Energy consumption [NOK]')
xlabel('Time [h]')
title('Cost function, J')
legend('J_{Optimal solution}','J_{Switching rule}')
hold off

%Saving in tikz
%matlab2tikz('minimumPointsDistance',{0.01},'optimal_solution_J.tikz')

%Plot of broad peak, optimal solution and solution from switching rule 

figure(3)
subplot(5,1,1)
hold on
stairs(t_p_broadpeak,P_broadpeak*60*60,'b')
ylabel('Energy price [kr/kWh]')
xlabel('Time [h]')
title('Energy price, P')
hold off

subplot(5,1,2)
hold on
plot(t_broadpeak,Tf_broadpeak,'b')
plot(t_broadpeak_s,Tf_broadpeak_s,'r')
plot(t_broadpeak,ones(length(t_broadpeak),1)*Tfmax-273,'k --')
plot(t_broadpeak,ones(length(t_broadpeak),1)*Tfmin-273,'k --')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tf')
legend('Optimal solution','Switching rule','T_{max}','T_{min}')
hold off

subplot(5,1,3)
hold on
plot(t_broadpeak,Qhf_broadpeak,'b')
plot(t_broadpeak_s,Qhf_broadpeak_s,'r')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhf') 
hold off

subplot(5,1,4)
hold on
plot(t_broadpeak,Tr_broadpeak,'b')
plot(t_broadpeak_s,Tr_broadpeak_s,'r')
plot(t_broadpeak,ones(length(t_broadpeak),1)*Trmax-273,'k --')
plot(t_broadpeak,ones(length(t_broadpeak),1)*Trmin-273,'k --')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tr')  
hold off

subplot(5,1,5)
hold on
plot(t_broadpeak,Qhr_broadpeak,'b')
plot(t_broadpeak_s,Qhr_broadpeak_s,'r')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhr')
hold off

%Saving in tikz
%matlab2tikz('minimumPointsDistance',{0.01},'broad_peak_TfTruP.tikz')

figure(4)
hold on
plot(t_broadpeak,J_broadpeak,'b')
plot(t_broadpeak_s,J_broadpeak_s,'r')
ylabel('Energy consumption [NOK]')
xlabel('Time [h]')
title('Cost function, J')
legend('J_{Optimal solution}','J_{switching rule}')
 hold off

 %Saving in tikz
%matlab2tikz('minimumPointsDistance',{0.01},'broad_peak_J.tikz')

%Plot of broad peak 2: no time to heat up the floor
%Optimal solution and solution from switching rule 

figure(5)
subplot(5,1,1)
stairs(t_p_broadpeak2,P_broadpeak2*60*60,'b')
ylabel('Energy price [kr/kWh]')
xlabel('Time [h]')
title('Energy price, P')

subplot(5,1,2)
hold on
plot(t_broadpeak2,Tf_broadpeak2,'b')
%plot(t_broadpeak2_s,Tf_broadpeak2_s,'r')
plot(t_broadpeak2,ones(length(t_broadpeak2),1)*Tfmax-273,'k --')
plot(t_broadpeak2,ones(length(t_broadpeak2),1)*Tfmin-273,'k --')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tf')
legend('Optimal solution','Switching rule','T_{max}','T_{min}')
hold off

subplot(5,1,3)
hold on
plot(t_broadpeak2,Qhf_broadpeak2,'b')
%plot(t_broadpeak2_s,Qhf_broadpeak2_s,'r')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhf') 
hold off

subplot(5,1,4)
hold on
plot(t_broadpeak2,Tr_broadpeak2,'b')
%plot(t_broadpeak2_s,Tr_broadpeak2_s,'r')
plot(t_broadpeak2,ones(length(t_broadpeak2),1)*Trmax-273,'k --')
plot(t_broadpeak2,ones(length(t_broadpeak2),1)*Trmin-273,'k --')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tr')  
hold off

subplot(5,1,5)
hold on
plot(t_broadpeak2,Qhr_broadpeak2,'b')
%plot(t_broadpeak2_s,Qhr_broadpeak2_s,'r')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhr')
hold off

%Saving in tikz
%matlab2tikz('minimumPointsDistance',{0.01},'case_3_TrTruP.tikz')

figure(6)
hold on
plot(t_broadpeak2,J_broadpeak2,'b')
%plot(t_broadpeak2_s,J_broadpeak2_s,'r')
ylabel('Energy consumption [NOK]')
xlabel('Time [h]')
title('Cost function, J')
legend('J_{Optimal solution}','J_{Switching rule}')
hold off

%Saving in tikz
%matlab2tikz('minimumPointsDistance',{0.01},'case_3_J.tikz')