%Optimal Temperature Control of rooms for Minimum Energy Cost

%Master thesis 2013

close all 
clear all
clc

% Single shooting optimization
% 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;

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

%Disturbances
To = 273-0; % [K]

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

%% Energy price

%Varying price
ts = 600; % 10min
tend = 54000; % 15h
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(61)); % 6 - 10
t3 = length(t(62):ts:t(end)); % 10 - end

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

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

par.P = P; 

%% Solving equations

%Constant parameters
% u = [qhf qhr]; %Constant u
d = [To min]; %Constant d


%Time varying parameters
% Input, u
qhf_min = 0; % [kW]
qhf_max = 2.5; % [kW]
qhr_min = 0; % [kW]
qhr_max = 2; %4.8; % [kW]

t1 = length(t(1):ts:t(25)); % 0:10min:4h
t2 = length(t(26):ts:t(37)); % 4:10min:6
t3 = length(t(38):ts:t(end)); % 6:10min:end

qhf1 = [ones(t1,1)*qhf_min; ones(t2,1)*qhf_max; ones(t3,1)*qhf_min]; % 15h
qhr1 = [ones(t1,1)*qhr_min; ones(t2,1)*qhr_max; ones(t3,1)*qhr_min]; % 15h

tintp = [0:ts:tend]'; 

u = [qhf1 ; qhr1];

[t x] = ode15s(@(t,x)equations(t,x,u,d,tintp,par),t,init);

%% Optimization
 
% Upper and lower constraints for the states
Tr_min = 293;  % [K]
Tr_max = 300;  % [K]
Tf_min = 292; % [K]
Tf_max = 300; % [K]

par1.Tr_min = Tr_min;
par1.Tr_max = Tr_max;
par1.Tf_min = Tf_min;
par1.Tf_max = Tf_max;

%Upper and lower bounds for the input
qhf_upper = 2.5; % [kW] 
qhr_upper = 2; %4.8; % [kW]
qhf_lower = 0; % [kW] 
qhr_lower = 0; % [kW] 

LB = [ones(length(t),1)*qhf_lower ; ones(length(t),1)*qhr_lower]; 
UB = [ones(length(t),1)*qhf_upper ; ones(length(t),1)*qhr_upper]; 

%Initial values
% load heatfile.mat  %Get the values from model_master_simulink
% qhf_g = interp1(time_sim,qhf_sim,t);
% qhr_g = interp1(time_sim,qhr_sim,t);
% u0 = [qhf_g ; qhr_g];
% u0 = [qhf_sim ; qhr_sim];
u0 = [2.5*ones(length(t),1); 2*ones(length(t),1)]; % [kW] qhf qhr

options = optimset('Display','Iter','Algorithm','active-set',...
    'TolCon',4*10^(-2),'TolX',0.01); 
Z = fmincon(@(u)costfun(u,t,init,d,tintp,par),u0,[],[],[],[],LB,UB,...
    @(u)nonlinconst(u,d,tintp,par1,par,t,init),options);


qhf_opt = Z(1:length(Z)/2,1); 
qhr_opt = Z(length(Z)/2+1:end,1);
[t x] = ode15s(@(t,x)equations(t,x,Z,d,tintp,par),t,init);

% Cost function
J = P.*(Z(1:length(Z)/2,1)+Z(length(Z)/2+1:end,1));

