% Script for finding the measurement combination matrices (H) for the 
% and exact local methods with the dynamic process model. 
% The setpoints of c=Hy are also calcualted

close all
clc
warning off

%%
open('CSTRandColumnAConnected')             % Open the Simulink model
sim('CSTRandColumnAConnected')              % Run the Simulink model

% Storing the optimal values of the measurements
y0=[T.signals.values(end,1);T.signals.values(end,6);T.signals.values(end,10);
    T.signals.values(end,14);T.signals.values(end,18);T.signals.values(end,22);
    ws.signals.values(end,10);ws.signals.values(end,11);ws.signals.values(end,2);
    ws.signals.values(end,3);ws.signals.values(end,5);ws.signals.values(end,12)];

%% Null space method in Simulink
% Disturbing the process by a 1% change in F0
set_param('CSTRandColumnAConnected/F0_step','before','0','after','460/60*0.01')
dF0=460/60*0.01;                            % Magnitude of input change

sim('CSTRandColumnAConnected')              % Run the Simulink model

% Storing the new values of the measurements
y_stepF0=y_n.signals.values(end,:)';

% Removing the disturbance in F0
set_param('CSTRandColumnAConnected/F0_step','before','0','after','460/60*0.01*0')

% Calculating the optimal sensitivity matrix
F=(y_stepF0-y0)./dF0;                       % Optimal sensitivity matrix
Hnsb=null(F');
Hnsb=Hnsb';                                 % Measurement combination matrix
Hns=Hnsb(3,:);                              % Choosing one of the rows
cs_ns=Hns*y0;                               % Set-point of c=H*y
% csns=335.8686                             % Calculated value of setpoint cs=H*y0

%% Exact local method in Simulink
% Disturbing the process by a 1% change in L
set_param('CSTRandColumnAConnected/LT - step','before','0','after','778/60*0.01')
sim('CSTRandColumnAConnected')              % Run the Simulink model

stepL=778/60*0.01;                          % Magnitude of input change
y_stepL=y_n.signals.values(end,:)';         % New measurement values

GyL=(y_stepL-y0)/stepL;                     % Finding dy/dL

% Removing the change in L
set_param('CSTRandColumnAConnected/LT - step','before','0','after','0')

% Disturbing the process by a 1% change in F
set_param('CSTRandColumnAConnected/Fs_step','before','0','after','958/60*0.01')
sim('CSTRandColumnAConnected')              % Run the Simulink model

stepF=958/60*0.01;                          % Magnitude of input change
y_stepF=y_n.signals.values(end,:)';         % New measurement values
GyF=(y_stepF-y0)/stepF;                     % Finding dy/dF

% Removing the change in F
set_param('CSTRandColumnAConnected/Fs_step','before','0','after','0')

% Disturbing the process by a 1% change in V
set_param('CSTRandColumnAConnected/VB_step','before','0','after','1276.1/60*0.01')

sim('CSTRandColumnAConnected')              % Run the Simulink model

stepV=1276.1/60*0.01;                       % Magnitude of input change
y_stepV=y_n.signals.values(end,:)';         % New measurement values
GyV=(y_stepV-y0)/stepV;                     % Finding dy/dF

% Removing the change in V
set_param('CSTRandColumnAConnected/VB_step','before','0','after','0')

% Computing the steady-state gain matrix Gy
Gy=[GyL GyF GyV];
% 2% measurement noise in flows
p=0.02;                                     
% Magnitude of noise in temperatures was chosen as 0.01
% Cimputing the scaling matrixes for noise and disturbance (Wn and Wd)
Wn=diag([0.01; 0.01; 0.01; 0.01; 0.01; 0.01; p*y0(7:12)]);  % 
F0=460/60;
% Magnitude of disturbance in F0 was chosen as 1%
dd=0.01*F0*ones(12,1);
Wd=diag(dd);

Y=F'*Wd*Wn;                             % Y matrix
Helb=(inv(Y*Y')*Gy)';                   % Measurement combination matrix
% Helb=((Y.\(Y'))*Gy)';
HelS=Helb(1,:);                         % Choosing one of the rows
cs_el=HelS*y0;                          % Set-point of c=H*y
%cs_el=-1.0e4;                          % Calculated value of setpoint cs=H*y0

