clc
clear
close all

N=10; % cells
M = 100; % simulations

F_c = 10; % cold fluid mass flow, kg/s


F_c_max = 10; % maximum cold fluid mass flow, kg/s
F_h_max = 3;  % maximum hot fluid mass flow, kg/s 
z_h_max = 1;
F_h = 0; % hot fluid mass flow,  kg/s

z_c = F_c/F_c_max;  % cold valve 
z_h = F_h./F_h_max;  % hot valve

rho_c = 980 ;        % cold fluid density, kg/m3
rho_h = 890;         % hot fluid density,  kg/m3

Tin_h = 70+273.15 ; % hot fluid inlet temperature, K
U = 150*1e-3;       % heat transfer coefficient,   kW/m2/K
A = 90;             % heat transfer area, m2
cp_c = 1.5;         % specific heat for cold fluid, kJ/kg/K
cp_h = 1.2;         % specific heat for hot fluid,  kJ/kg/K

Tin_c = linspace(22.7,26.3,M)+273.15;
p = [rho_c; rho_h; Tin_h; U; A; cp_c; cp_h;N;F_c_max;F_h_max]; % parameters

Tc0 = [314.1411  310.3426  307.0824  304.2841  301.8823  299.8208  298.0514  296.5326  295.2291  294.1102]';
Th0 = [335.2388  328.4495  322.6229  317.6224  313.3307  309.6473  306.4859  303.7724  301.4434  299.4442]';
x0 = [Tc0; Th0];
ThSP = x0(end);


Kc_c = 0;
tauI_c = 1;
Kc_h = 0.0803;
tauI_h = 86;


for i=1:M    
u = [z_c; z_h]; % input
d = Tin_c(i);
sim('TestControllerV2')
z_h_save(i) = u2(end);
end

%% Save Data
save('output.mat','z_h_save','Tin_c','N')
clear
load('output.mat')



%% Plot saved data
figure
plot(z_h_save,Tin_c-273.15,'k','Linewidth',1)
hold on
plot([1, 1], [min(Tin_c), max(Tin_c)]-273.15, 'k:','Linewidth',1)
ylabel('Cold inlet temperature,\textsuperscript{o}C','Interpreter','latex')
xlabel('Hot fluid valve opening','Interpreter','latex')
ax = gca;
ax.FontSize = 16;
ax.TickLabelInterpreter = 'latex';

% Add annotation
startB1 = [0.7535 0.7169]; % point of starting annotation
endB1 = [0.7534 0.5349]; % ending annotation
B1 = annotation('textarrow',startB1,endB1,'String','Bottleneck 1 - cooling water valve','Interpreter','latex','FontSize',16); % bottlenneck 1 - cooling water

startB2 = [0.8 0.5];
endB2 = [1 0.5];
B2 = annotation('textarrow', startB2,endB2,'String','Bottleneck 2 - hot fluid valve valve','Interpreter','latex','FontSize',16); % bottlenneck 2 - hot fluid

posIR = [0.5 0.6 0.3 0.3]; % position of infeasible region text box in normalized coordinates
posFR = [0.2 0.1 0.3 0.3]; % position of feasible region text box in normalized coordinates
IR = annotation('textbox',posIR,'String','Infeasible Region','FitBoxToText','on','Interpreter','latex','FontSize',16); % ingeasible region
FR = annotation('textbox',posFR,'String','Feasible Region','FitBoxToText','on','Interpreter','latex','FontSize',16);   % feasible region


