function runblending()

% Run the blending process
% Martha Grover, July 6, 2017

% Here is the playlist for the accompanying screencasts:
% https://www.youtube.com/playlist?list=PL4xAk5aclnUhb0tM6nypIATyxPRk0fB3L

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Dynamic Simulation of the Nonlinear System %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Following on the first screencast (Blending Process: Dynamic Simulation),
% here you can simulate the nonlinear coupled set of differential equations.
% As in the second screencast (Blending Process: Steady States), 
% observe here how the mass fraction in the tank 
% approaches the steady state value of xbar as the system approaches long 
% time.  The volume does not have a unique steady state.  Since filling a 
% tank is a integrating process, it may increase or decrease without bound, 
% based on the difference between the inlet and outlet mass flow rates.

clc
clear all
close all

% Parameters
global rho 
rho = 1000;   % Density [kg/m^3] 
tf = 3600;    % Final simulation time [s]

% Initial conditions
V0 = 10;        % Initial volume in the tank [m^3]
x0 = 0.025;     % Initial mass fraction of Species A in the tank (unitless)
% Note: these are currently set at their nominal steady-state values...
% But you can try changing them to other values.
z0 = [V0; x0];  % Initial condition vector for the system

% Nominal inputs: 
% x1: Mass fraction of Species A in Stream 1 (unitless)
% w1: Mass flow rate in Stream 1 [kg/s]
% x2: Mass fraction of Species A in Stream 2 (unitless)
% w2: Mass flow rate in Stream 2 [kg/s]
% w:  Outlet mass flow [kg/s]

% Here we will consider x1 and x2 to remain as constants
global x1 x2
x1 = 0.1;
x2 = 0;
% The flow rates will be the variables in u that may change over time
w1 = 1;         % w1: Mass flow rate in Stream 1 [kg/s]
w2 = 3;         % w2: Mass flow rate in Stream 2 [kg/s]
w = w1 + w2;    % w:  Outlet mass flow [kg/s]
u = [w1; w2; w];    % Assemble these together in a column vector u

% Simulate the blending process
[t,z] = ode45(@blending, [0 tf],z0,[],u);

% Plot the blending process under steady operation
figure(1)
subplot(2,1,1)
plot(t,z(:,1))
grid on
xlabel('time t [s]')
ylabel('V [m^3]')
title('Blending process under constant inputs')
subplot(2,1,2)
plot(t,z(:,2))
grid on
xlabel('time t [s]')
ylabel('x')

% Calculate the steady-state value for x 
xbar = (w1*x1 + w2*x2)/(w1+w2);
% Check to see if this matches the long-term value of x on the plot
% Since V does not have a unique steady state, specify a desired operating value
% Here it is set equal to the initial volume
Vbar = V0;  % volume [m^3]
% Define ubar as the set of inputs to achieve this steady state
ubar = u;

% Next simulate the system for the higher value of w1 = 1.1 (0.1 higher)
dw1 = 0.1;
u = ubar + [dw1; 0; 0];
[tstep,zstep] = ode45(@blending, [0 tf],z0,[],u);
% Plot the blending process under a step change in w1
figure(2)
subplot(2,1,1)
plot(tstep,zstep(:,1))
grid on
xlabel('time t [s]')
ylabel('V [m^3]')
title('Blending process under a step input to w_1')
subplot(2,1,2)
plot(tstep,zstep(:,2))
grid on
xlabel('time t [s]')
ylabel('x')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% Linear Approximate Model %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Formulate the linear approximation as a state-space system using
% deviation variables (See screencast on Blending Process: Deviation
% Variables)

% Dynamic state x: V' and x'

% Inputs u: w1', w2', w'

% Measurements y: w1', V', x'

% Define the matrices A, B, C, and D for the linearized system
% dxdt = A*x + B*u' 
% y = C*x + D*u'
% See the screencasts on Linearization and State Space

% Dynamic equation: dx/dt = Ax + Bu
global A B
A = [0 0; 0 1/(rho*Vbar)*(-w1-w2)];   
B = [1/rho 1/rho -1/rho; 1/(rho*Vbar)*(x1-xbar) 1/(rho*Vbar)*(x2-xbar) 0];
% Measurement equation: y = Cx + Du
C = [0 0; 1 0; 0 1];            % Measuring the states V and x directly
D = [1 0 0; 0 0 0; 0 0 0];      % Measuring the input w directly

% Define input in deviation variables u'
up = u - ubar;
% Define the deviation in initial state
% Here we will continue to use z for the state vector, to avoid confusion 
% with the mass fraction of Species A which is called x
zbar = [Vbar; xbar];
zp0 = z0 - zbar;   % Note: We set z0 = zbar earlier, so this should be zero
% Simulate the linearized model of the blending process
[tp,zp] = ode45(@blending_linear,[0 tf],zp0,[],up);

% Plot the comparison between the nonlinear system and the linearized
% system, for a step change of 0.1 in w1
figure(3)
subplot(2,1,1)
plot(tstep,zstep(:,1),tp,zp(:,1)+Vbar)   
grid on
legend('Nonlinear','Linearized')
xlabel('time t [s]')
ylabel('volume V [m^3]')
title('Comparison of nonlinear and linear models for a step change in w_1')
subplot(2,1,2)
plot(tstep,zstep(:,2),tp,zp(:,2)+xbar)   
grid on
legend('Nonlinear','Linearized')
xlabel('time t [s]')
ylabel('mass fraction x')

% In this example, the linear model agrees pretty well with the nonlinear 
% model.  The equation for volume V is linear and does not depend on the 
% mass fraction x, so it agrees perfectly with the nonlinear model.
% The mass fraction prediction is not perfect, but because the step change
% is "small" (10%) the agreement is still pretty good.

% Steady-state gain of the system
K = C*inv(-A)*B+D;
% Note: It is not possible to calculate this since A is not invertible
% This makes sense because there is no unique steady-state value for
% volume V

% Note that we did not actually use C and D in this example.  We can
% evaluate the measurements y' at each time point of our ode solution:
yp = zeros(length(ubar),length(tp));
for i = 1:length(tp)
   yp(:,i) = C*zp(i,:)' + D*up;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Define the Nonlinear and Linear Models %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dzdt = blending(t,z,u)
% Dynamic model of the blending process
% Martha Grover, July 6, 2017
% State variables passed in as vector z:
V = z(1);   % Volume [m3]
x = z(2);   % Mass fraction of Species A in the tank (unitless)
% Mass fractions in Streams 1 and 2: x1 and x2
global x1 x2
% Inputs passed in as vector u:
w1 = u(1);  % Mass flow rate in Stream 1 [kg/s]
w2 = u(2);  % Mass flow rate in Stream 2 [kg/s]
w = u(3);   % Here we set the outlet flow equal to the total inlet flow

% Density rho
global rho

% Dynamic equations
dVdt = 1/rho*(w1+w2-w);                     % Change in volume
dxdt = 1/(rho*V)*(w1*(x1-x)+w2*(x2-x));     % Change in mass fraction

% Assemble the two time derivatives back into a vector to pass back
dzdt = [dVdt; dxdt];

end

function dxdt = blending_linear(t,x,u)
% State space version of the linearized model for the blending process,
% using deviation variables
% dx/dt = A*x + B*u
% The measurement equation of y = C*x + D*u is not included here, since
% it is an algebraic equation that does not need to be solved by ode45.
% Martha Grover, July 6, 2017

% Dynamic state x
Vp = x(1);      % V' = deviation in volume
xp = x(2);      % x' = deviation in mass fraction

% Input vector u 
w1p = u(1);     % w1' = deviation in mass flow rate in inlet Stream 1
w2p = u(2);     % w2' = deviation in mass flow rate in inlet Stream 2
wp = u(3);      % w' = deviation in mass fraction in outlet stream

% Use the values in A and B defined in the main script
global A B

% dx/dt = Ax + Bu
dxdt = A*x + B*u;

end






