%% In this function you:
% - Define the differential equations for your model
% - Linearize and get your model in the the state-space form.
% - Obtain the transfer function matrices G(s) and Gd(s)
%
% To obtain the transfer functions G and Gd from the state-space form:
% Suppose you have the linearized system x_dot=Ax+Bu+Ed, y=C*x; Eq.(1)
% where x=delta states, u= delta MV, d= delta DV and y= delta CV
% From this we want to compute the transfer function matrices G(s) and Gd(s)
% 
% We apply Laplace transform to Eq (1) to obtain: 
% G(s)=C*(I*s-A)^-1*B and Gd(s)=C*(I*s-A)^-1*E  
% You should know how to derive this by hand! possible exam question 
%
% You can also consider that u inclues MVs and DVs. In this case the model
% would be x_dot=Ax+Bd*u, where Bd is a nx-by-(nMV+nDV) matrix. 
% In this case you would get only one transfer function matrix G(s). 
% Both forms are equivalent, as long you know what your are doing.
% 
%%

function [G,Gd,A,B,C,D,E]=computeTransferFunctions


%%Definition of differential equations

% Define/create symbolic variables for your MV, DV, CV
syms h q1 q2 cA cAf k

% Tip: if your equations have additional variables that depend on
% your MVs, DVs or CVs and you want to define these additional variables, 
% so that you can write your differential equations in a simpler way, 
% you can declare these "aid" variables as symbolic variables too. For
% examples: valve_flow, enthalpy_in, etc.

%Define vectors of DV(d), MV(u), states(x) and CV(y)
d = [cAf;k];           %deviation variables; vertical vector (if more than 1 DV)
u = [q1;q2];                        %inputs; vertical vector (if more than 1 MV)
x = [h;cA];                  %states; vertical vector (if more than one states)
y = [h;cA];                     % outputs; vertical vector (if more than one CV)
% y can be different than the states x, but we recommend to select y=x as 
% in part 3 you may decide to control all your states


%Parameters that will be used in the differential equations
Area = 4;  %m^2 - Note that we are not calling it A because below A will be  
          % matrix A from x_dot=Ax+Bu+Ed, y=C*x. In other words, be careful
                                              % when naming your parameters.

%Here you can write additional algebraic equations for the variables that 
%you will use in the differential equations. 

%Differential equations
f(1,1) = 1/Area*(q1-q2);
f(1,2) = (cAf - cA)*(q1/(Area*h)) - k*cA^2;
%For more than 1 state: f(1,1)=...; f(1,2)=...;
%Alternatively, write f directly as an horizontal vector

%% Linearization
% Define the Jacobian - for the state space form x_dot=Ax+Bu+Ed
A=jacobian(f,x); 
B=jacobian(f,u);
E=jacobian(f,d);

% Define the Jacobian from states to output: y=Cx+Du
C=jacobian(y,x);                       % Select the outputs from the states 
D=jacobian(y,u);                % Matrix from input to output. Usually  D=0

%Define the nominal point for x,y, u, d 
h = 1;
cA = 0.05;

q1 = 1;
q2 = 1;

cAf = 1;
k = 95;

%Replace the symbolic variables by their corresponding nominal value
A=double(subs(A));
B=double(subs(B));
E=double(subs(E));

C=double(subs(C));
D=double(subs(D));

%Laplace variable
s=tf('s');

%Identity matrix
n=length(A); I=eye(n);

%Obtaining transfer function matrix - y = G*u+Gd*d
% you must know how to do this by hand! possible exam question.
G=C*inv(s*I-A)*B+D;                 % tranfer function from input to output
Gd=C*inv(s*I-A)*E;            % tranfer function from disturbance to output

%Simplifying equation - minimum realization (cancels common roots in numerator and denominator)
G=minreal(G);
Gd=minreal(Gd);

%You may change the display format.
G=set(zpk(G), 'DisplayFormat','time constant');                   
Gd=set(zpk(Gd), 'DisplayFormat','time constant'); 
end

