%Optimal Temperature Control of rooms for Minimum Energy Cost

%Master thesis 2013

close all 
clear all
clc

%Optimal control using PI

%Function to find the tuning parameters of the system
%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 = 24*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;
%var = [t; P]';
%var.signals.dimensions = dim;

%% Inputs and initial conditions

%Constant inputs
%qhr = 2; % No heat on: 0 %max: 2 [kW] 
%qhf = 2.5; % No heat on: 0 %max: 2.5 [kW] 

%Disturbance
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 = 293; % [K]
J = 0; % [NOK]
init = [Tf Tr J];

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

%Upper and lower constraints
Tfmax = 298; % [K]
Tfmin = 293; % [K] 293
Trmax = 298; % [K]
Trmin = 292; % [K] 293

%% Simulation without PI controller

sim('simulink_mastermdl')

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

%Simulation time
t = time/(60*60); %[h] 

% Plotting
figure(1)
Subplot(4,1,1)
plot(t,Tf,'b')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tf') 

subplot(4,1,2)
plot(t,Qhf,'b')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhf')

Subplot(4,1,3)
plot(t,Tr,'r')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tr') 

subplot(4,1,4)
plot(t,Qhr,'r')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhr')

figure(2)
plot(t,J,'g')
ylabel('Cost [kr]')
xlabel('Time [h]')
title('Cost function, J')

%% Simulation with PI controller

%Simulations to find the tuning parameters

sim('simulink_master_PI')

%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)

%Plot to find the tuning parameters
figure(3)
subplot(2,1,1)
plot(t,Tf,'b')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tf')

subplot(2,1,2)
plot(t,Qhf,'b')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhf') 

%Saving in tikz
%matlab2tikz('open_loop_Tf_plot.tikz')

figure(4)
subplot(2,1,1)
plot(t,Tr,'b')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tr')

subplot(2,1,2)
plot(t,Qhr,'b')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhr') 

%Saving in tikz
%matlab2tikz('open_loop_Tr_plot.tikz')

%% Simulation of control structure
sim('simulink_master_PI_control_structure')

%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)

%Plot from control structure

figure(5)
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)
plot(t,Tf,'b')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tf')

subplot(5,1,3)
plot(t,Qhf,'b')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhf') 

subplot(5,1,4)
plot(t,Tr,'b')
ylabel('Temperature [C]')
xlabel('Time [h]')
title('Tr')  

subplot(5,1,5)
plot(t,Qhr,'b')
ylabel('Heat [kW]')
xlabel('Time [h]')
title('qhr')

figure(6)
plot(t,J,'b')
ylabel('Energy consumption [NOK]')
xlabel('Time [h]')
title('Cost function, J')
