%Optimal Temperature Control of rooms for Minimum Energy Cost

%Master thesis 2013
%Siri Hofstad Trapnes

close all 
clear all
clc

% Function for finding the solution of optimization problem. 
% Written by Chriss Grimholt, 2013


%% Parameters
par.UA_fr = 0.350; %[kJ/sK = kW/K]
par.UA_ro = 0.007; %[kJ/sK = kW/K]

par.mCp_f = 4000;  %[kJ/K]
par.mCp_r = 70;    %[kJ/K]
par.m_r   = 70;    %[kg]

% Simulation parameters
opt.sim_time = 20; %[h]

% Input Bounds
opt.lb.u.q_f = 0;   %[kW]
opt.ub.u.q_f = 2.5; %[kW]

opt.lb.u.q_r = 0;   %[kW]
opt.ub.u.q_r = 2;   %[kW]

% State Bounds
opt.lb.x.T_r = 19;  %[C]
opt.ub.x.T_r = 25;  %[C]

opt.lb.x.T_f = 20;  %[C]
opt.ub.x.T_f = 25;  %[C]

% Initial States
init.T_f  = 20;    %[C]
init.T_r  = 19;    %[C]  
init.J    = 0;     %[Kr] 
init.vector = [init.T_f init.T_r init.J 0 0]';

% Disturbances
dist.price.value = [0.6 1 0.6];  %[kr/kWh]
dist.price.time  = [0 6 8];     %[h]

dist.m_in.value = [0.06];   %[kg] 
dist.m_in.time  = [0];      %[h] 

%Constant outdoor temperture
dist.T_o.value = 0;    %[C] 
dist.T_o.time  = 0;    %[h] 
%Varying outdoor temperature
% 1) step before the peak period:
% dist.T_o.value = [0 -9];  %[C]
% dist.T_o.time = [0 0];    %[h]
% 2) step after the peak period
% dist.T_o.value = [0 -9];  %[C]
% dist.T_o.time = [0 9];    %[h]

%% Optimization

%Storing in room and floor
% Initial guess:
%  qhf qhr
% u = [2 2];   
% u = fminsearch(@(u) costfun(u,dist,par,init,opt), u)

% %Storing in floor
% Initial guess:
u = 2.91; %Time
u = fminsearch(@(u) costfun([u 15],dist,par,init,opt), u)
% %where 15 is a number after t(p_high)

%% Simulation
%u = [u 10] %Optimal value
u = [(6-6.01) 10] % Optimal value switching rule
option=odeset('AbsTol',1e-10,'RelTol',1e-3);
[t, x] = ode23(@(t,x) model(t,x,u,dist,par,opt),[0 opt.sim_time],init.vector(1:3),option);
cost = x(end,end) %The cost


%% Plotting

% Back calculating u
T_o = zeros(length(t),1);
m_in = zeros(length(t),1);
for i = 1:length(t)
   
    T_o_tmp = dist.T_o.value(dist.T_o.time <= t(i));
    T_o(i) = T_o_tmp(end);

    m_in_tmp = dist.m_in.value(dist.m_in.time <= t(i));
    m_in(i) = m_in_tmp(end);
end
u_f  = zeros(length(t),1);
u_r  = zeros(length(t),1);

t_int = find((t < u(1))>0);
u_f(t_int)  = max(par.UA_fr*(x(t_int,1)-x(t_int,2)),0);
t_int = find((t < u(2))>0);
u_r(t_int) = max(-par.UA_fr*(x(t_int,1)-x(t_int,2)) + (par.UA_ro + m_in(t_int)./par.m_r*par.mCp_r).*(x(t_int,2)-T_o(t_int)),0);

t_int = find((t > u(1)) .* (t <= dist.price.time(2))>0);
u_f(t_int) = opt.ub.u.q_f; %#ok<FNDSB>
t_int = find((t > u(2)) .* (t <= dist.price.time(2))>0);
u_r(t_int) = opt.ub.u.q_r;  %#ok<FNDSB>

t_int = find((t > u(1)) .* (t <= dist.price.time(2) .* (x(:,1) >= opt.ub.x.T_f))>0);
u_f(t_int)  = max(par.UA_fr*(x(t_int,1)-x(t_int,2)),0);
t_int = find((t > u(2)) .* (t <= dist.price.time(2) .* (x(:,2) >= opt.ub.x.T_r))>0);
u_r(t_int) = max(-par.UA_fr*(x(t_int,1)-x(t_int,2)) + (par.UA_ro + m_in(t_int)./par.m_r*par.mCp_r).*(x(t_int,2)-T_o(t_int)),0);


t_int = find(((t >dist.price.time(2)) .* (x(:,1) <= opt.lb.x.T_f))>0);
u_f(t_int) = max(par.UA_fr*(x(t_int,1)-x(t_int,2)),0);

t_int = find(((t > dist.price.time(2)) .* (x(:,2) <= opt.lb.x.T_r))>0);
u_r(t_int) =max( -par.UA_fr*(x(t_int,1)-x(t_int,2)) + (par.UA_ro + m_in(t_int)./par.m_r*par.mCp_r).*(x(t_int,2)-T_o(t_int)),0);


% Result from constant disturbances

figure(1)
subplot(411)
plot(t,x(:,1:2))
xlabel('Time,[h]')
ylabel('Temp,[C]')
ylim([18,25])
subplot(412)
plot(t,[u_f,u_r]);
legend('Floor','Room')
xlabel('Time, [h]')
ylabel('Power,[kW]')
subplot(413)
plot(t,x(:,3))
xlabel('Time, [h]')
ylabel('Cost, [kr.]')
subplot(414)
stairs([dist.price.time, opt.sim_time],[dist.price.value, dist.price.value(end)])
xlabel('Time, [h]')
ylabel('Price,[kr.]')

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

%% Results from varying disturbances

%Saving variables for disturbances before the peak
% savefile = 'disturbance9_bp.mat'
% x_9bp = x;
% t_9bp = t;
% To_9bp = T_o;
% u_r_9bp = u_r;
% u_f_9bp = u_f;
% save(savefile,'x_9bp','To_9bp','u_r_9bp','u_f_9bp','t_9bp')
% savefile = 'disturbance5_bp.mat'
% x_5bp = x;
% t_5bp = t;
% To_5bp = T_o;
% u_r_5bp = u_r;
% u_f_5bp = u_f;
% save(savefile,'x_5bp','To_5bp','u_r_5bp','u_f_5bp','t_5bp')
% savefile = 'disturbance_minus5_bp.mat'
% x_m5bp = x;
% t_m5bp = t;
% To_m5bp = T_o;
% u_r_m5bp = u_r;
% u_f_m5bp = u_f;
% save(savefile,'x_m5bp','To_m5bp','u_r_m5bp','u_f_m5bp','t_m5bp')
% savefile = 'disturbance_minus9_bp.mat'
% x_m9bp = x;
% t_m9bp = t;
% To_m9bp = T_o;
% u_r_m9bp = u_r;
% u_f_m9bp = u_f;
% save(savefile,'x_m9bp','To_m9bp','u_r_m9bp','u_f_m9bp','t_m9bp')

% %Saving variables for disturbances after the peak period
% savefile = 'disturbance9_afterpeak.mat'
% x_9ap = x;
% t_9ap = t;
% To_9ap = T_o;
% u_r_9ap = u_r;
% u_f_9ap = u_f;
% save(savefile,'x_9ap','To_9ap','u_r_9ap','u_f_9ap','t_9ap')
% savefile = 'disturbance5_afterpeak.mat'
% x_5ap = x;
% t_5ap = t;
% To_5ap = T_o;
% u_r_5ap = u_r;
% u_f_5ap = u_f;
% save(savefile,'x_5ap','To_5ap','u_r_5ap','u_f_5ap','t_5ap')
% savefile = 'disturbance_minus5_afterpeak.mat'
% x_m5ap = x;
% t_m5ap = t;
% To_m5ap = T_o;
% u_r_m5ap = u_r;
% u_f_m5ap = u_f;
% save(savefile,'x_m5ap','To_m5ap','u_r_m5ap','u_f_m5ap','t_m5ap')
% savefile = 'disturbance_minus9_afterpeak.mat'
% x_m9ap = x;
% t_m9ap = t;
% To_m9ap = T_o;
% u_r_m9ap = u_r;
% u_f_m9ap = u_f;
% save(savefile,'x_m9ap','To_m9ap','u_r_m9ap','u_f_m9ap','t_m9ap')

%Saving variables for disturbances before the peak. Switching rules
% savefile = 'disturbance9_bps.mat'
% x_9bps = x;
% t_9bps = t;
% To_9bps = T_o;
% u_r_9bps = u_r;
% u_f_9bps = u_f;
% save(savefile,'x_9bps','To_9bps','u_r_9bps','u_f_9bps','t_9bps')
% savefile = 'disturbance5_bps.mat'
% x_5bps = x;
% t_5bps = t;
% To_5bps = T_o;
% u_r_5bps = u_r;
% u_f_5bps = u_f;
% save(savefile,'x_5bps','To_5bps','u_r_5bps','u_f_5bps','t_5bps')
% savefile = 'disturbance_minus5_bps.mat'
% x_m5bps = x;
% t_m5bps = t;
% To_m5bps = T_o;
% u_r_m5bps = u_r;
% u_f_m5bps = u_f;
% save(savefile,'x_m5bps','To_m5bps','u_r_m5bps','u_f_m5bps','t_m5bps')
% savefile = 'disturbance_minus9_bps.mat'
% x_m9bps = x;
% t_m9bps = t;
% To_m9bps = T_o;
% u_r_m9bps = u_r;
% u_f_m9bps = u_f;
% save(savefile,'x_m9bps','To_m9bps','u_r_m9bps','u_f_m9bps','t_m9bps')

% %Result from varying disturbances before peak 
% load disturbance9_bp.mat
% load disturbance5_bp.mat
% load disturbance_minus5_bp.mat
% load disturbance_minus9_bp.mat

% figure(2)
% subplot(411)
% hold on
% plot(t,x(:,1),'b')
% plot(t_9bp,x_9bp(:,1),'r')
% plot(t_5bp,x_5bp(:,1),'g')
% plot(t_m5bp,x_m5bp(:,1),'y')
% plot(t_m9bp,x_m9bp(:,1),'m')
% xlabel('Time,[h]')
% ylabel('Temp,[C]')
% title('T_F')
% ylim([18,25])
% hold off
% 
% subplot(412)
% hold on
% plot(t,x(:,2),'b')
% plot(t_9bp,x_9bp(:,2),'r')
% plot(t_5bp,x_5bp(:,2),'g')
% plot(t_m5bp,x_m5bp(:,2),'y')
% plot(t_m9bp,x_m9bp(:,2),'m')
% xlabel('Time,[h]')
% ylabel('Temp,[C]')
% legend('constant To','plus9','plus5','minus5','minus9')
% ylim([18,25])
% title('T_R')
% hold off
% 
% subplot(413)
% hold on
% plot(t,u_f,'b');
% plot(t_9bp,u_f_9bp,'r');
% plot(t_5bp,u_f_5bp,'g');
% plot(t_m5bp,u_f_m5bp,'y');
% plot(t_m9bp,u_f_m9bp,'m');
% title('q_{H/F}')
% xlabel('Time, [h]')
% ylabel('Power,[kW]')
% hold off
% 
% subplot(414)
% hold on
% plot(t,u_r,'b');
% plot(t_9bp,u_r_9bp,'r');
% plot(t_5bp,u_r_5bp,'g');
% plot(t_m5bp,u_r_m5bp,'y');
% plot(t_m9bp,u_r_m9bp,'m');
% title('q_{H/R}')
% xlabel('Time, [h]')
% ylabel('Power,[kW]')
% hold off
% % 
% %Saving in tikz
%  matlab2tikz('case1_dist_at0_Tu.tikz')
% 
% figure(3)
% % subplot(311)
% % stairs([dist.price.time, opt.sim_time],[dist.price.value, dist.price.value(end)])
% % xlabel('Time, [h]')
% % ylabel('Price,[kr.]')
% 
% subplot(211)
% hold on
% plot(t,T_o,'b')
% plot(t_9bp,To_9bp,'r')
% plot(t_5bp,To_5bp,'g')
% plot(t_m5bp,To_m5bp,'y')
% plot(t_m9bp,To_m9bp,'m')
% legend('constant To','plus9','plus5','minus5','minus9')
% xlabel('Time [h]')
% ylabel('Temp [C]')
% title('Outdoor temperature')
% hold off
% 
% subplot(212)
% hold on
% plot(t,x(:,3),'b')
% plot(t_9bp,x_9bp(:,3),'r')
% plot(t_5bp,x_5bp(:,3),'g')
% plot(t_m5bp,x_m5bp(:,3),'y')
% plot(t_m9bp,x_m9bp(:,3),'m')
% %legend('constant To','plus9','plus5','minus5','minus9')
% xlabel('Time [h]')
% ylabel('Cost [kr]')
% hold off
% % 
% % % Saving in tikz
%   matlab2tikz('case1_dist_at0_PJTo.tikz')

%Result from varying disturbances after peak 
% load disturbance9_afterpeak.mat
% load disturbance5_afterpeak.mat
% load disturbance_minus5_afterpeak.mat
% load disturbance_minus9_afterpeak.mat
% 
% figure(2)
% subplot(411)
% hold on
% plot(t,x(:,1),'b')
% plot(t_9ap,x_9ap(:,1),'r')
% plot(t_5ap,x_5ap(:,1),'g')
% plot(t_m5ap,x_m5ap(:,1),'y')
% plot(t_m9ap,x_m9ap(:,1),'m')
% xlabel('Time,[h]')
% ylabel('Temp,[C]')
% title('T_F')
% ylim([18,25])
% hold off
% 
% subplot(412)
% hold on
% plot(t,x(:,2),'b')
% plot(t_9ap,x_9ap(:,2),'r')
% plot(t_5ap,x_5ap(:,2),'g')
% plot(t_m5ap,x_m5ap(:,2),'y')
% plot(t_m9ap,x_m9ap(:,2),'m')
% xlabel('Time,[h]')
% ylabel('Temp,[C]')
% legend('constant To','plus9','plus5','minus5','minus9')
% ylim([18,25])
% title('T_R')
% hold off
% 
% subplot(413)
% hold on
% plot(t,u_f,'b');
% plot(t_9ap,u_f_9ap,'r');
% plot(t_5ap,u_f_5ap,'g');
% plot(t_m5ap,u_f_m5ap,'y');
% plot(t_m9ap,u_f_m9ap,'m');
% title('q_{H/F}')
% xlabel('Time, [h]')
% ylabel('Power,[kW]')
% hold off
% 
% subplot(414)
% hold on
% plot(t,u_r,'b');
% plot(t_9ap,u_r_9ap,'r');
% plot(t_5ap,u_r_5ap,'g');
% plot(t_m5ap,u_r_m5ap,'y');
% plot(t_m9ap,u_r_m9ap,'m');
% title('q_{H/R}')
% xlabel('Time, [h]')
% ylabel('Power,[kW]')
% hold off
% 
% %Saving in tikz
% matlab2tikz('case1_dist_at9_Tu.tikz')
% 
% figure(3)
% % subplot(311)
% % stairs([dist.price.time, opt.sim_time],[dist.price.value, dist.price.value(end)])
% % xlabel('Time, [h]')
% % ylabel('Price,[kr.]')
% 
% subplot(211)
% hold on
% plot(t,T_o,'b')
% plot(t_9ap,To_9ap,'r')
% plot(t_5ap,To_5ap,'g')
% plot(t_m5ap,To_m5ap,'y')
% plot(t_m9ap,To_m9ap,'m')
% %legend('constant To','plus9','plus5','minus5','minus9')
% xlabel('Time [h]')
% ylabel('Temp. [C]')
% title('Outdoor temperature')
% hold off
% 
% subplot(212)
% hold on
% plot(t,x(:,3),'b')
% plot(t_9ap,x_9ap(:,3),'r')
% plot(t_5ap,x_5ap(:,3),'g')
% plot(t_m5ap,x_m5ap(:,3),'y')
% plot(t_m9ap,x_m9ap(:,3),'m')
% legend('constant To','plus9','plus5','minus5','minus9')
% xlabel('Time, [h]')
% ylabel('Cost, [kr.]')
% hold off
% 
% %Saving in tikz
%  matlab2tikz('case1_dist_at9_PJTo.tikz')

% Result from switching rule. Disturbance before peak period
% load disturbance9_bps.mat
% load disturbance5_bps.mat
% load disturbance_minus5_bps.mat
% load disturbance_minus9_bps.mat %Infeasible solution
% 
% figure(2)
% subplot(411)
% hold on
% plot(t,x(:,1),'b')
% plot(t_9bps,x_9bps(:,1),'r')
% plot(t_5bps,x_5bps(:,1),'g')
% plot(t_m5bps,x_m5bps(:,1),'y')
% %plot(t_m9bps,x_m9bps(:,1),'m') %Infeasible solution
% xlabel('Time,[h]')
% ylabel('Temp,[C]')
% title('T_F')
% ylim([18,25])
% hold off
% 
% subplot(412)
% hold on
% plot(t,x(:,2),'b')
% plot(t_9bps,x_9bps(:,2),'r')
% plot(t_5bps,x_5bps(:,2),'g')
% plot(t_m5bps,x_m5bps(:,2),'y')
% %plot(t_m9bps,x_m9bps(:,2),'m') %Infeasible solution
% xlabel('Time,[h]')
% ylabel('Temp,[C]')
% legend('constant To','plus9','plus5','minus5','minus9')
% ylim([18,25])
% title('T_R')
% hold off
% 
% subplot(413)
% hold on
% plot(t,u_f,'b');
% plot(t_9bps,u_f_9bps,'r');
% plot(t_5bps,u_f_5bps,'g');
% plot(t_m5bps,u_f_m5bps,'y');
% %plot(t_m9bps,u_f_m9bps,'m'); %Infeasible solution
% title('q_{H/F}')
% xlabel('Time, [h]')
% ylabel('Power,[kW]')
% hold off
% 
% subplot(414)
% hold on
% plot(t,u_r,'b');
% plot(t_9bps,u_r_9bps,'r');
% plot(t_5bps,u_r_5bps,'g');
% plot(t_m5bps,u_r_m5bps,'y');
% %plot(t_m9bps,u_r_m9bps,'m'); %Infeasible solution
% title('q_{H/R}')
% xlabel('Time, [h]')
% ylabel('Power,[kW]')
% hold off
% 
% % %Saving in tikz
% % matlab2tikz('case1_dist_at0_Tu_switch.tikz')
% 
% figure(3)
% subplot(311)
% stairs([dist.price.time, opt.sim_time],[dist.price.value, dist.price.value(end)])
% xlabel('Time, [h]')
% ylabel('Price,[kr.]')
% 
% subplot(312)
% hold on
% plot(t,T_o,'b')
% plot(t_9bps,To_9bps,'r')
% plot(t_5bps,To_5bps,'g')
% plot(t_m5bps,To_m5bps,'y')
% %plot(t_m9bps,To_m9bps,'m') %Infeasible solution
% legend('constant To','plus9','plus5','minus5','minus9')
% xlabel('Time [h]')
% ylabel('Temp. [C]')
% legend('Outdoor temperature')
% hold off
% 
% subplot(313)
% hold on
% plot(t,x(:,3),'b')
% plot(t_9bps,x_9bps(:,3),'r')
% plot(t_5bps,x_5bps(:,3),'g')
% plot(t_m5bps,x_m5bps(:,3),'y')
% %plot(t_m9bps,x_m9bps(:,3),'m') %Infeasible solution
% legend('constant To','plus9','plus5','minus5','minus9')
% xlabel('Time, [h]')
% ylabel('Cost, [kr.]')
% hold off
% 
% % % %Saving in tikz
% % %  matlab2tikz('case1_dist_at0_PJTo_switch.tikz')