% Script for calculating gains 
clear all; close all; clc;
global coldata dyndata hydraulic termodata UU alfa tend NI i
%% 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 %%
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
%Change in disturbances can be mase here to obtain Gyd
dyndata.x_nom = x_nom;        %initial composition and holdup states
% disturbance in feed F, h=1e-5
h = 1e-5;
%dd1= h;
%u_nom(6)= u_nom(6)+dd1;
%disturbance in composition A in feed
%dd2=h;
%u_nom(7)=u_nom(7)+dd2;
%disturbance in composition B in feed
%dd3=h;
%u_nom(8)=u_nom(8)+dd3;

%disturbance in boilup
%dd4=h;
%u_nom(2)=u_nom(2)+dd4;

%disturbance in vapor split RV
%dd5=h;
%u_nom(4)=u_nom(4)+dd5;


dyndata.u_nom=u_nom;                  %we have made initial input global
dyndata.temp_nom=temp_nom;
%% Run this file to get initial values with LV Configuartion, uncomment
% line number 
NI=5; % Number of integerators
options = odeset('MaxStep',10);
tend=20000;
[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)';
Plot_Profiles
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
%Temperature as measurement Y
[sys Temp]= Petlyuk_Nlmodel(t,xf,u_nom);
%Yss = Temp;                          % Measurements, Temperatures on all
%stages at steady state
%plot_temp_profile                   % Temperature Profile in the column

%% Calculate steady state gain matrix for three remaining D.O.F. (L,S,Rl) GY 
load Y0.mat              % same as nominal temperature
% For the disturbance in feed flow rate F
%Run the simulation and  store the value of Temp as YF 
%YF= Temp;
%dYdF = (YF-Y0)/dd1;
load YF.mat
load dYdF.mat

%For disturbance in composition of A in feed, zA
%Run the simulation and store the value of Temp as YzA
%YzA = Temp;
%dYdzA = (YzA-Y0)/dd2;
load YzA.mat
load dYdzA.mat

%For disturbance in composition of B in feed, zB
%Run the simulation and store the value of Temp as YzB
%YzB = Temp;
%dYdzB = (YzB-Y0)/dd3;
load YzB.mat
load dYdzB.mat

%For disturbance in the boilup V
%Run the simulation and store the value of Temp as YV
%YV = Temp;
%dYdV = (YV-Y0)/dd4;
load YV.mat
load dYdV.mat

%For disturbance in the vapor split ratio RV
%Run the simulation and store the value of Temp as YV
%YRV = Temp;
%dYdRV = (YRV-Y0)/dd5;
load YRV.mat
load dYdRV.mat

% Calculate the steady state gain model for diturbances
GYd = [dYdF dYdzA dYdzB dYdV dYdRV];

% Load steady-state disturbance gain
%load GYd
