% Script for calculating gain matrix GY (w.r.t. input) 

clear all; close all; clc;
global coldata dyndata hydraulic termodata UU alfa tend NI i stages_Temperature stages_Temperature_name Temp_nom
%% COLUMN DATA %%

% Number of stages in 6 sections of Petlyuk 
coldata.Nt = [12 12 12 12 12 12];
% Initial holdup stages
coldata.M0 = [1 1 1 1 1 1 1 1 1 1 1 1];
%Boiling points of ethanol, propanol and n-butanol
coldata.Tboil = [78.37 97 117.4];
termodata.P = 1.0;
hydraulic.k1 = 1/0.063;
%% THERMODYNAMIC DATA %%  (Realtive volatilities alfa)
alfa=[4.215374389	2.102194688	1] %ethanol propanol butanol

%% load nominal x u
load Nominal_final.mat        % The nominal value obtained from optimization for ...

%minimizing the sum of impurities in the product streams (D, S and R)
UU=[];                        %
i=0;%to initiate Petlyuk_openloop make change in CS by i+1
indmat_Yes_S1
% Make the nominal values of x, u and Temp global
dyndata.x_nom = x_nom; 
%u_nom(6) =0.9*u_nom(6);
%u_nom(2) = 1.1*u_nom(2);
dyndata.u_nom=u_nom;          
%dyndata.temp_nom=temp_nom;
%get stage indices for prefrac and main column
coldata.stages_pre=[im_x1(1:coldata.Nt(1)) im_x2(1:coldata.Nt(2))];
coldata.stages_main=[im_x21(1:coldata.Nt(3)) im_x22(1:coldata.Nt(4)) im_x23(1:coldata.Nt(5)) im_x24(1:coldata.Nt(6))];
%define position for temperature controller
stages_Temperature_name={'L','S1','RL'};
stages_Temperature=[im_x22(2) im_x23(12) im_x1(10)];
%calculate nominal temperature
Temp_nom=get_Temperature(stages_Temperature,1,dyndata.x_nom);
%% 

NI=5; % Number of integerators
options = odeset('MaxStep',10);
tend=5000;
[t,x]=ode15s(@Simulation_Petlyuk,[0 tend],[dyndata.x_nom;zeros(NI,1)],options);
%% xf states to plot profiles.
xf=x(end,1:225)';
u_nom=UU(end,2:12)';
dxdt = norm(Petlyuk_NLmodel(0, xf, u_nom));
% Print the values of compositions in D (distillate), S (side stream), R (residue)...
%and D1 and R1, (top and bottom product from the prefractionator both  end)
dis_composition
[sys Temp] = Petlyuk_NLmodel(tend, xf,u_nom);

% Optimal temperatures (measurement)also saved as temp_nom in
% Nominal_final.mat
load Y0.mat

%% Calculate steady state gain matrix for three remaining D.O.F. (L,S,Rl) GY 
%load Y0.mat
% Go to Petlyuk_openloop and give the step change in input L,uncomment line
% 53, 54
%  and comment line 57, 58 dyndata.u_nom(1)+1e-5
%Run the simulation  and  store the value of Temp as YL  
%YL= Temp;
%delL= 1e-5;
%dYdL = (YL-Y0)/delL;
load YL.mat
load dYdL.mat

%Give step change in S1, dyndata.u_nom(5)+1e-5
%Run the simulation and store the value of Temp as YS
%YS=Temp;
%delS= 1e-5;
%dYdS = (YS-Y0)/delS;
load YS.mat
load dYdS.mat
%Go to file Petlyuk_openloop and comment line number and uncomment line 
%Give step change in Rl, dyndata.u_nom(3)+1e-5
%Run the simulation and store the value of Temp as YRl
%YRl=Temp;
%delRl=1e-5;
%dYdRl = (YRl-Y0)/delRl;
load YRl.mat
load dYdRl.mat

%Steady state gain matrix GY
GY = [dYdL dYdS dYdRl];
