% Main script for simulating the room temperature control using SRC
% September 2018

%There are 4 available MVs
% 1: Refrigerated water
% 2: Cooling water
% 3: Hot water
% 4: Electric heating

%Here we only use the calculated Kc and alphas from the system of
%equations- 

clc
clear
close all

%Load results from simulation
%a  % slopes
%Kc % overall controller gain
%du % delta u for each MV
%load('SRCvalues_Kc_alphaV3'); %design values from solving the system of equations that mantains Kuy constant
load('SRCvalues_Kc_alpha_delay'); %design values from solving the system of equations that mantains Kuy constant
%% Simulation parameters
tend=250;           % Simulation time

%% Nominal values
Tr0=   18;           % Room temperature
Tamb0= 18;           % Ambient temperature
Tsp0=  18;           % Set point temperature


v10=0;              % v1 opening
v20=0;  
v30=0;
v40=0;

u0=0;               % Controller output; %from the SRC solution, CW crosses 0 at 43.7500

%% Disturbances

% Set point
Tsp1= 5;            % First set point
Tsp2=8;            % Second change of set point
Tsp3= -Tsp1-Tsp2-2;            % Third change of set point
Tsp4= -9; 

tTsp1= 10;          % First set point, time
tTsp2= tTsp1+50;    % Second change of set point, time
tTsp3= tTsp2+50;    % Third change of set point, time
tTsp4= tTsp3+50; 

% Ambient temperature
Tamb1=0;
Tamb2=0;
Tamb3=0;
Tamb4=0;

tTamb1=tTsp1;
tTamb2=tTsp2;
tTamb3=tTsp3;
tTamb4=tTsp4;


%% Upper and lower limits

% Max and minimum values for valve openings. This is NOT 0-1 because here
% we are working in deviation variables. 
% Still not sure if we actually need this. For the time being, -0.5 to 0.5
% v1max=0.5; 
% v1min=-v1max;
% v2max=0.5;
% v2min=-v2max;
% v3max=0.5;
% v3min=-v3max;
% v4max=0.5;
% v4min=-v4max;

v1max=vmax; 
v1min=0;
v2max=vmax;
v2min=0;
v3max=vmax;
v3min=0;
v4max=vmax;
v4min=0;

% For each MV, the following limits are the split points in the SRC
% diagram. In the Simulink file we are working in deviation variables. 
split0=[c100(1), c100(2), c0(2), c100(3), c100(4)];
split=split0-split0(3)*ones(1,length(split0));


RWlow= split(1);
RWup= split(2);
CWlow= split(2);
CWup= split(3);
HWlow= split(3);
HWup= split(4);
EHlow=split(4);
EHup= split(5);

%% Correction considering offsets and scaling in SRC
%Scaling of signals for v
% This has to be done because in the system of equations we are solving, we
% do not assure that this condition is fulfilled. 

Ks1= vmax/du(1);
Ks2=vmax/du(2);
Ks3=vmax/du(3);
Ks4=vmax/du(4);
Ksu=1; %1/100;

% We should consider offsets
off1=-split(2);
off2=0;
off3=0;
off4=-split(4);



%% Transfer functions from valves to room temperature

%Refrigerated water
K_rw=k(1);   %gain must be negative! more RW, lower temperature
tau_rw=tau(1);
theta_rw=theta(1);

%Cooling water
K_cw=k(2);    %gain must be negative! more CW, lower temperature
tau_cw=tau(2);
theta_cw=theta(2);

%Heating water
K_hw=k(3);
tau_hw=tau(3);
theta_hw=theta(3);

%Electric heating
K_eh=k(4);
tau_eh=tau(4);
theta_eh=theta(4);

%Ambient temperature (Disturbance)
K_amb=1;
tau_amb=15;
theta_amb=6;

taup=[tau_rw,tau_cw,tau_hw,tau_eh];
thetap=[theta_rw, theta_cw, theta_hw, theta_eh];

%% Controller tunning for SRC
% Kc is calculated from the system of equations. 
%Here we only have to calculate tauI

taup_tune=mean(taup);
theta=mean(thetap); 
tauc=mean(tauC);

tauI_sim=min(taup_tune,4*(tauc+theta));
tauI_sim=max(tauI);
%Kc=Kc;
%tauI=4; %vpa(tauI);
%% Corrections for valve openings 

% Gain correction  because of the way that the SRC is calculated, the sign
% of the gain has to be given somewhere.
K1=sign(Kcp(1));
K2=sign(Kcp(2));
K3=sign(Kcp(3));
K4=sign(Kcp(4));



% %% Simulations with changes in set-points
% sim('SRC_Room_simpleV3');%,'simulationTime', [0 tend])  
%   
% %Figures
% disp(' Set points')
% disp('  IAE no lead-lag')
% IAE(end)
% 
% 
% set(0,'DefaultTextFontName','Times',...
% 'DefaultTextFontSize',30,...
% 'DefaultAxesFontName','Times',...
% 'DefaultAxesFontSize',30,...
% 'DefaultLineLineWidth',2.5,...
% 'DefaultLineMarkerSize',12)
% 
% figure  
% subplot(2,1,1)
% plot(t,Troom,'r'); %'Color',[ 138/255,43/255,226/255])
% hold on
% plot(t,Tsp,'k--')
% plot(t,Tamb,'Color',[34/255,139/255,34/255])
% %grid on
% %grid minor
% ax = gca;
% ax.FontSize = 16;
% %xlabel('Time,s')
% ylim([0 40])
% ylabel('Temperature, T (°C)','Interpreter','Latex')
% legend('T','T^{ref}','T^{amb}','Location','Northeast','Orientation','horizontal')
% hold off
% %legend boxoff
% 
% subplot(2,1,2)
% plot(t,v1,'Color',[25/255    25/255    112/255])
% hold on
% %plot(t,v2,'Color',[0    0.4470    0.7410])
% plot(t,v2,'Color',[0.3010    0.7450    0.9330])
% 
% % plot(t,v1,'Color',[0.3010    0.7450    0.9330])
% % hold on
% % plot(t,v2,'Color',[0    0.4470    0.7410])
% plot(t,v3,'Color', [0.9290    0.6940    0.1250])
% plot(t,v4,'r')
% %plot(t,u,'k'); %'LineStyle',':','Color',[0,128/255,0])
% %grid on
% %grid minor
% ax = gca;
% ax.FontSize = 16;
% %title('Valve openings to process - without lead-lag','Interpreter','Latex')
% xlabel('Time (min)')
% ylim([-0.05 1.05])
% ylabel('Valve opening, $u_i$  (-)','Interpreter','Latex')
% %ylabel('$v_i$, u','Interpreter','Latex')
% legend('u_{AC}','u_{CW}','u_{HW}','u_{EH}', 'Location','Northeast','Orientation','horizontal')
% %ylim([-10 100])
% %xlim([0 tend])
% hold off
% %legend boxoff

%

%% Simulations with changes in ambient temperature

% Tsp1= 0;            % First set point
% Tsp2=0;            % Second change of set point
% Tsp3=0;            % Third change of set point
% Tsp4=0; 
% 
% % Ambient temperature
% Tamb1=2;
% Tamb2=10;
% Tamb3=-13;
% Tamb4=-15;

tTamb1=10;
tTamb2=80;
tTamb3=140;
tTamb4=180;
tTamb5=280;
tTamb6=350;

% Ambient temperature
Tamb1=2; %0.2;
Tamb2=9; %0.10;
Tamb3=-5; %-0.12;
Tamb4=-25; %-0.15;
Tamb5=18;
Tamb6=5;
tend=400;

sim('SRC_Room_simpleV4');%,'simulationTime', [0 tend])                
 
disp('  Ambient temperature')
disp('  IAE no lead-lag')
IAE(end)



figure  %Only without lead lag
subplot(2,1,1)
plot(t,Troom,'r');%'Color',[ 138/255,43/255,226/255])
hold on
%plot(t_ll,Troom_ll,'k'); %'Color',[ 138/255,43/255,226/255])
plot(t,Tsp,'k--')
plot(t,Tamb,'Color',[34/255,139/255,34/255])
%grid on
%grid minor
ax = gca;
ax.FontSize = 16;
%xlabel('Time,s')
%ylim([0 40])
ylabel('Temperature, T (°C)','Interpreter','Latex')
legend('T','T^{ref}','T^{amb}','Location','Northeast','Orientation','horizontal')
%title('Changes in $T^{amb}$','Interpreter','Latex')
hold off
%legend boxoff

subplot(2,1,2)
plot(t,v1,'Color',[25/255    25/255    112/255])
hold on
%plot(t,v2,'Color',[0    0.4470    0.7410])
plot(t,v2,'Color',[0.3010    0.7450    0.9330])
plot(t,v3,'Color', [0.9290    0.6940    0.1250])
plot(t,v4,'r')
%plot(t,u,'k'); %'LineStyle',':','Color',[0,128/255,0])
%grid on
%grid minor
ax = gca;
ax.FontSize = 16;
%title('Valve openings to process - without lead-lag','Interpreter','Latex')
xlabel('Time (min)')
ylim([-0.05 1.05])
%ylabel('$v_i$, u','Interpreter','Latex')
%legend('v_{AC}','v_{CW}','v_{HW}','v_{EH}','u', 'Location','Southeast','Orientation','horizontal')
ylabel('Valve opening, $u_i$ (-)','Interpreter','Latex')
legend('u_{AC}','u_{CW}','u_{HW}','u_{EH}', 'Location','Northeast','Orientation','horizontal')
%ylim([-10 100])
%xlim([0 tend])
hold off
%legend boxoff

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