Contents

 clear; clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution for ChE 231 MATLAB Assignment 3 %
%              Spring 2017                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Given Parameters

Define given parameters

Y_XS=0.4;  %g/g
alpha=2.5;  %g/g
beta=0.15;  %/hr
mu_m=0.8;   %/hr
Pm=52;      %g/L
Km=1.7;     %g/L
Ki=20;      %g/L
D=0.2;     %/hr
Sf=10;      %g/L
p=[Y_XS alpha beta mu_m Pm Km Ki D Sf];

% % Plot differential equations at various concentrations
% % in order to visually find range around each steady state
% X_vec=0:0.1:4; S_vec=0:0.1:2; P_vec=0:1:5;
%
% for k=1:length(P_vec)
%     for j=1:length(S_vec)
%         for i=1:length(X_vec)
%             X=X_vec(i);
%             S=S_vec(j);
%             P=P_vec(k);
%             F=bioreactor([X,S,P],p);
%             dX(i,j)=F(1);
%             dS(i,j)=F(2);
%             dP(i,j)=F(3);
%         end
%     end
%     figure(k)
%     mesh(X_vec,S_vec,dX','k')
%     hold on
%     mesh(X_vec,S_vec,dS')
%     mesh(X_vec,S_vec,dP')
%     legend('dX','dS','dP')
%     xlabel('X'); ylabel('S'); zlabel('differential');
%     clear dX dS dP
% end

Question 1

display('----------------------------------------------')
display('Question 1:')

% solves washout steady state
x0=[0 100 0];
options=optimoptions('fsolve','display','off');
F_SS1=fsolve(@(x) bioreactor(x,p),x0,options);
display(F_SS1)

% solves washout steady state
x0=[1 1 1];
options=optimoptions('fsolve','display','off');
F_SS2=fsolve(@(x) bioreactor(x,p),x0,options);
display(F_SS2)
display('----------------------------------------------')
----------------------------------------------
Question 1:

F_SS1 =

    3.6669    0.8327   11.9175


F_SS2 =

    3.6669    0.8327   11.9175

----------------------------------------------

Question 2

display('Question 2:')
display('Linearization dx/dt=Ax, where x=[X S P]'':')
syms X S P
f=bioreactor([X S P],p);
A=jacobian(f,[X S P]);
X=F_SS2(1); S=F_SS2(2); P=F_SS2(3);
A=subs(A);
A=double(A);
display(A);
display('Eigen values of A:')
lambda=eig(A);
display(lambda)
display('----------------------------------------------')
Question 2:
Linearization dx/dt=Ax, where x=[X S P]':

A =

   -0.0000    0.5713   -0.0183
   -0.5000   -1.6283    0.0457
    0.6500    1.4283   -0.2457

Eigen values of A:

lambda =

   -1.4719
   -0.2022
   -0.2000

----------------------------------------------

Question 3

display('Question 3:')
display('(See MATLAB code for implementation of ode45 function)')
tspan=[0 10];
[t,sol]=ode45(@(t,x) bioreactor(x,p),tspan,F_SS2);

display('----------------------------------------------')
Question 3:
(See MATLAB code for implementation of ode45 function)
----------------------------------------------

Question 4

display('Question 4:')
display('(See figure 1)')
figure(1)
plot(t,sol)
legend('X','S','P'); xlabel('time (hr)'); ylabel('Concentration (g/L)');
title('Steady state dynamic response');
display('The dynamic responses are straight lines because the system')
display('is initially at steady state. Differentials equal zero.')
display('----------------------------------------------')
Question 4:
(See figure 1)
The dynamic responses are straight lines because the system
is initially at steady state. Differentials equal zero.
----------------------------------------------

Question 5

display('Question 4:')
D=0.15;     %/hr
p=[Y_XS alpha beta mu_m Pm Km Ki D Sf];
[t,sol]=ode45(@(t,x) bioreactor(x,p),tspan,F_SS2);
display('(See figure 2)')
figure(2)
plot(t,sol)
legend('X','S','P'); xlabel('time (hr)'); ylabel('Concentration (g/L)');
title('Dynamic response with D=0.15 hr^{-1}');

D=0.25;     %/hr
p=[Y_XS alpha beta mu_m Pm Km Ki D Sf];
[t,sol]=ode45(@(t,x) bioreactor(x,p),tspan,F_SS2);
display('(See figure 3)')
figure(3)
plot(t,sol)
legend('X','S','P'); xlabel('time (hr)'); ylabel('Concentration (g/L)');
title('Dynamic response with D=0.25 hr^{-1}');
set(gca,'FontSize',12); set(gca,'FontName','Times');
display('Biomass concentration increases with lower dilution rate,')
display('and decreases with higher dilution rate.')
display('----------------------------------------------')
Question 4:
(See figure 2)
(See figure 3)
Biomass concentration increases with lower dilution rate,
and decreases with higher dilution rate.
----------------------------------------------

Question 6

display('Question 4:')
% First find the new jacobian matrix using new diltion rate
D=0.15;     %/hr
p=[Y_XS alpha beta mu_m Pm Km Ki D Sf];
syms X S P
f=bioreactor([X S P],p);
A=jacobian(f,[X S P]);
X=F_SS2(1); S=F_SS2(2); P=F_SS2(3);
A=subs(A);
A=double(A);
% Solve system linearized using jacobian matrix
[t,sol]=ode45(@(t,x) bioreactor_linear(x,A),tspan,F_SS2);
display('(See figure 4)')
figure(4)
plot(t,sol)
legend('X','S','P'); xlabel('time (hr)'); ylabel('Concentration (g/L)');
title('Linearized dynamic response with D=0.15 hr^{-1}');

% First find the new jacobian matrix using new diltion rate
D=0.25;     %/hr
p=[Y_XS alpha beta mu_m Pm Km Ki D Sf];
syms X S P
f=bioreactor([X S P],p);
A=jacobian(f,[X S P]);
X=F_SS2(1); S=F_SS2(2); P=F_SS2(3);
A=subs(A);
A=double(A);
% Solve system linearized using jacobian matrix
[t,sol]=ode45(@(t,x) bioreactor_linear(x,A),tspan,F_SS2);
display('(See figure 5)')
figure(5)
plot(t,sol);
legend('X','S','P'); xlabel('time (hr)'); ylabel('Concentration (g/L)');
title('Linear dynamic response with D=0.25 hr^{-1}');
display('Similar response in short time, however, non-realistic')
display('results with substrate concentration going negative.')
display('----------------------------------------------')
Question 4:
(See figure 4)
(See figure 5)
Similar response in short time, however, non-realistic
results with substrate concentration going negative.
----------------------------------------------