% Valve with actuator.
% May 2019

% This script contains parameters for the steady state and dynamic part of
% the valve. The ss part is the valve equation that calculates flow as 
% a function of opening.  The dynamic part is a 2nd order ODE. 
% The script also calculates the tuning for the PI controller that controls  
% the MeOH mass flow, using the input to the valve actuator as MV. 

%There are two parts of the valve: the actuator and the actual body/stem. 
%The actuator provides the dynamics using a system of 2 ODEs. 

%Functions that use what is generated here:
%The model for the actuator is in: actuator.m
%The model for the valve (flow as function of opening) is in valve.m

%The stem provides the actual flow using the valve characteristic (static
%equation). 
% Possible valve characteristics: "valve_type": 
    %1=> equal percentage
    %2=> square root
    %3=> linear
  
%Functions that use what is generated here:
%The model for the actuator is in: actuator.m 
%The model for the stem is in: valve.m

%In the formaldehyde case study, we will use the same model for the 
% actuator for both the valve that feeds MeOH and the blower that feeds 
% the air.  For the blower, the actual MV is the opening of the internal 
% guide vanes (IGV), and it is assumed that the dynamic behavior is similar
% to a valve. 
% The difference will be the model of the unit and the dynamic parameters 
% of the actuator.


%% Parameters for MeOH valve

%The labels on the right (comments), correspond to the names in the test
%scripts. I changed them to avoid confusions between MeOH and air flows and
%paramters (especially regarding the actuator). 

m_MeOH_sp_0=2940;               % kg/h %Flow_sp0
m_MeOH_max=4200;                %kg/h; %Flow_max
m_MeOH_min=0;                   %kg/h
valve_delta_p=0.5;              %Nominal pressure drop; eventually this could be a DV %delta_p
valve_Cv=m_MeOH_max/sqrt(valve_delta_p);

% Valve characteristic: "valve_type": 
    %1=> equal percentage
    %2=> square root
    %3=> linear
valve_type=3; %actually, the whole thing is tested only for a linear valve. Someone else in the future can test this with other valve characteristics.  

% Valve actuator parameters
valve_z0=0.7;               %initial opening (main state for the ODE)
valve_act_Ga=1;                   %actuator gain
valve_act_tau=1;                  %actuator time constant
valve_act_zeta=1.1;               %actuator damping factor

valve_w0=0;                       %initial condition of aid variable for ODE system. One equation is dS/dt=w, where S is the stem movement (z), and at initial conditions, it should not be moving. 

I_0_valve=valve_z0/valve_act_Ga;      %Initial controller output signal

valve_actuator_0=[valve_z0 valve_w0]; %Initial conditions for valve actuator
valve_actuator_par=[valve_act_Ga;valve_act_tau;valve_act_zeta];


%% Calculate transfer function and tune

%FOTD transfer function --- this is assuming overdamped system (zeta>1)
s=tf('s');
% p1=(-valve_act_zeta+sqrt(valve_act_zeta^2-1))/act_tau;
% p2=(-valve_act_zeta-sqrt(valve_act_zeta^2-1))/act_tau;
gvalve=valve_Cv/(valve_act_tau^2*s^2+2*valve_act_tau*valve_act_zeta*s+1);
gvalve = zpk(gvalve);
gvalve = set(gvalve,'DisplayFormat','time constant');
[z,p,k] = zpkdata(gvalve);
tau_valve=abs(1./p{1});

%half-rule approximation
tau1_valve=max(tau_valve);
tau2_valve=min(tau_valve);
K_valvefotd=valve_Cv;
theta_valvefotd=tau2_valve/2;           %2nd order tf has theta=0;
tau_valvefotd=tau1_valve+tau2_valve/2;   
gaprox_valve=K_valvefotd*exp(-theta_valvefotd*s)/(tau_valvefotd*s+1);
%step(gvalve,gaprox_valve)

%SIMC tuning
tauc_valve=theta_valvefotd*2;
Kc_valve=tau_valvefotd/(K_valvefotd*(tauc_valve+theta_valvefotd));
tauI_valve=min(tau_valvefotd,(4*(tauc_valve+tau_valvefotd)));
KI_valve=Kc_valve/tauI_valve;

