clear all
clc
%Interfacing to current UniSim file and identifing of unit operations
h=ActXServer('UnisimDesign.Application');
hyCase = h.Activedocument;
sol=hyCase.Solver;
f = hyCase.Flowsheet;
desorber = f.operation.Item('T-101');
desFS = desorber.ColumnFlowsheet;
desTray = desFS.Operations.Item('Main TS');
absorber = f.operation.Item('T-100');
absFS = absorber.ColumnFlowsheet;

% Base Case Data; FF:Feed Flowrate, TT:Tray Temperature 
BaseFF = f.MaterialStreams.Item('1').MolarFlowValue;
BaseCO2CompFF = f.MaterialStreams.Item('1').ComponentMolarFractionValue(1);
BaseCO2CompDes = f.MaterialStreams.Item('9').ComponentMolarFractionValue(1);
BaseDesTT = desTray.TemperatureValue;

% Make a change in compositions of a flowrate is not possible from Matlab so we define pure flow for creating disturbance in compositions
BaseN2 = f.MaterialStreams.Item('N2').MolarFlowValue;
BaseCO2 = f.MaterialStreams.Item('CO2').MolarFlowValue;
BaseO2 = f.MaterialStreams.Item('O2').MolarFlowValue;
BaseH2O = f.MaterialStreams.Item('H2O').MolarFlowValue;
BaseAmine = f.MaterialStreams.Item('Amine').MolarFlowValue;

% ***First Disturbance (Feed Flow rate:+5%); Fdist: Feed flowrateasdisturbance
sol.CanSolve=0;
Fdist = 1.05;
f.MaterialStreams.Item('N2').MolarFlowValue = BaseN2*Fdist;
f.MaterialStreams.Item('CO2').MolarFlowValue = BaseCO2*Fdist;
f.MaterialStreams.Item('O2').MolarFlowValue = BaseO2*Fdist;
f.MaterialStreams.Item('H2O').MolarFlowValue = BaseH2O*Fdist;
f.MaterialStreams.Item('Amine').MolarFlowValue = BaseAmine*Fdist;
invoke (desFS, 'Reset');
sol.CanSolve=1;
invoke (desFS, 'Run');

%******* Value of different parameters in the process after first disturbance******** 
% FDistFF: Feed flowrate with its new value as disturbance
% FDistCO2CompDes: CO2 composition in the bottom of stripper after first distrubance
% FDistDesTT: Desorber Trays Temperature after first distrubance
FDistFF = f.MaterialStreams.Item('1').MolarFlowValue;
FDistCO2CompDes = f.MaterialStreams.Item('9').ComponentMolarFractionValue(1);
FDistDesTT = desTray.TemperatureValue;

%*************Disturbance sensitivity *********************************
% DeltaCO2d1: difference of optimum CO2 composition in bottom of desorber after first distrurbance compare to base case
% DeltaTempd1: difference of optimum trays temperature of desorber after first distrurbance compare to base case
DeltaCO2d1 = FDistCO2CompDes - BaseCO2CompDes;
DeltaTempd1 = FDistDesTT - BaseDesTT;

%***** Second Disturbance (CO2 Composition in Feed:+5%)*******
sol.CanSolve=0;

% New Compositions
CO2Dif = BaseCO2*abs(Fdist-1);
f.MaterialStreams.Item('CO2').MolarFlowValue = BaseCO2*Fdist;
FlowSum =BaseN2 + BaseO2 + BaseH2O + BaseAmine;
% 5% increase in CO2 composition normalize on other compositions
f.MaterialStreams.Item('N2').MolarFlowValue = BaseN2 - CO2Dif * BaseN2 / FlowSum;
f.MaterialStreams.Item('O2').MolarFlowValue = BaseO2 - CO2Dif * BaseO2 / FlowSum;
f.MaterialStreams.Item('H2O').MolarFlowValue = BaseH2O - CO2Dif * BaseH2O / FlowSum;
f.MaterialStreams.Item('Amine').MolarFlowValue = BaseAmine - CO2Dif * BaseAmine / FlowSum;

invoke (desFS, 'Reset');
sol.CanSolve=1;
invoke (desFS, 'Run');

%******* Value of different parameters in the process after 2nd disturbance******** 
% CompDistCO2CompFF: CO2 compostion in feed with its new value as disturbance
% CompDistCO2CompDes: CO2 composition in the bottom of stripper after 2nd distrubance
% CompDistDesTT: Desorber Trays Temperature after 2nd distrubance
CompDistCO2CompFF = f.MaterialStreams.Item('1').ComponentMolarFractionValue(1);
CompDistCO2CompDes = f.MaterialStreams.Item('9').ComponentMolarFractionValue(1);
CompDistDesTT = desTray.TemperatureValue;

% Disturbance sensitivity
% DeltaCO2d2: difference of optimum CO2 composition in bottom of desorber after 2nd distrurbance compare to base case
% DeltaTempd2: difference of optimum trays temperature of desorber after 2nd distrurbance compare to base case
DeltaCO2d2 = CompDistCO2CompDes - BaseCO2CompDes;
DeltaTempd2 = CompDistDesTT - BaseDesTT;

%***** Third Disturbance (Change in Stripper pressure 10kPa)*******
% Back to Base Case
sol.CanSolve=0;
f.MaterialStreams.Item('N2').MolarFlowValue = BaseN2;
f.MaterialStreams.Item('CO2').MolarFlowValue = BaseCO2;
f.MaterialStreams.Item('O2').MolarFlowValue = BaseO2;
f.MaterialStreams.Item('H2O').MolarFlowValue = BaseH2O;
f.MaterialStreams.Item('Amine').MolarFlowValue = BaseAmine;
invoke (desFS, 'Reset');
sol.CanSolve=1;
invoke (desFS, 'Run');
% Change in pressure: Bottom from 200 kPa to 210 kPa and Top from 180 kPa to 190kPa
% PressDistCO2CompDes: CO2 composition in the bottom of stripper after 3rd distrubance
% PressDistDesTT: Desorber Trays Temperature after 3rd distrubance
f.MaterialStreams.Item('9').PressureValue=210;
f.MaterialStreams.Item('8').PressureValue=190;
invoke (desFS, 'Run');
PressDistCO2CompDes = f.MaterialStreams.Item('9').ComponentMolarFractionValue(1);
PressDistDesTT = desTray.TemperatureValue;

% Disturbance sensitivity
% DeltaCO2d3: difference of optimum CO2 composition in bottom of desorber after 3rd distrurbance compare to base case
% DeltaTempd3: difference of optimum trays temperature of desorber after 3rd distrurbance compare to base case
DeltaCO2d3 = PressDistCO2CompDes - BaseCO2CompDes;
DeltaTempd3 = PressDistDesTT - BaseDesTT;

% Overal Optimal Variation
% DeltaCO2Opt: overal optimal variation of CO2 composition in the bottom of desorber
% DeltaTempOpt: overal optimal variation of tray temperatures in desorber
DeltaCO2Opt = (DeltaCO2d1^2 + DeltaCO2d2^2 + DeltaCO2d3^2)^0.5;
DeltaTempOpt = (DeltaTempd1.^2 + DeltaTempd2.^2 + DeltaTempd3.^2).^0.5;

% Back to Base Case
sol.CanSolve=0;
f.MaterialStreams.Item('9').PressureValue=200;
f.MaterialStreams.Item('8').PressureValue=180;
invoke (desFS, 'Reset');
sol.CanSolve=1;
invoke (desFS, 'Run');

% ***************************Finding Gain Value****************************
% We want to find difference in CO2 composition in the bottom of desorber
% and difference in trays temperature of dersorber when small change in
% reboiler duty happens. While CO2 composition is sepc of the tower we
% change it and find change in reboiler duty and trays temperature profile
desSpec = desFS.ActiveSpecifications.Item('Comp Fraction');
BaseRebDuty = f.EnergyStreams.Item('Q-101').HeatFlowValue;
change=0.009;
DeltaCO2 = desSpec.GoalValue*change;
invoke (absFS, 'Reset');
desSpec.GoalValue = DeltaCO2+desSpec.GoalValue;
invoke (desFS, 'Reset');
invoke (absFS, 'Run');
invoke (desFS, 'Run');

RebDuty = f.EnergyStreams.Item('Q-101').HeatFlowValue;
DestrayTemp=desTray.TemperatureValue;
DeltaTemp = DestrayTemp - BaseDesTT;

% Difference in manipulate variable 
DeltaU = RebDuty - BaseRebDuty;

% Gain value for CO2 composition in the bottom of desorber
GCO2 = DeltaCO2/DeltaU;
% Gain value for tray temperatures of desorber
GTemp = DeltaTemp ./ DeltaU;

% IE:Impelementation Error
IECO2 = 0.001;
IETemp = 1;
SpanCO2 = DeltaCO2Opt + IECO2;
SpanTemp = DeltaTempOpt + IETemp;

% Scaled Gain
% SGCO2: Scaled Gain of CO2 composition; SGTemp: Scaled Gain of tray temperatures
SGCO2 = abs(GCO2/SpanCO2);
SGTemp = abs(GTemp ./ SpanTemp);





