function FedBatch4
clear, clc, format short g, format compact
Sin=90;
for k=1:20
    Sin=Sin+5;
    Sin_list(k)=Sin;
    PR(k)=-FedBatch3(Sin);
    disp([ ' Sin = ' num2str(Sin) ' Product = ' num2str(PR(k) )])
end
plot(Sin_list, PR)
xlabel(' Feed Glucose Concentration (g/L) ');
ylabel(' Production in 10 cycles (gr) '); 
pause
Sinopt= fminbnd(@FedBatch3,Sin_list(1), Sin_list(end)) ;
disp([ 'Soptim =  ' num2str(Sinopt) ' Product (gr) = ' num2str(FedBatch3(Sinopt))]);
function Product= FedBatch3 (Sin)
tspan = [0 1.]; % Range for the independent variable 
y0 = [24.; 0; 0; 0.8; 0]; % Initial values for the dependent variables 
tg=[]; Nx=[]; Ns=[];Np=[]; V=[]; PR=[];
[t,y]=ode45(@ODEfunIni,tspan,y0,[],Sin);
%disp([' time = '  int2str(t(end) ) '  Nx = '  num2str(y(end,1)) ' Ns = ' num2str(y(end,2)) ' Np = ' num2str(y(end,3)) ' V = ' num2str(y(end,4)) ' Pr= ' num2str(y(end,5))]); 
tspan=[t(end) t(end)+5];
y0=y(end,:);
tg=[tg; t]; Nx=[Nx; y(:,1)]; Ns=[Ns; y(:,2)]; Np=[Np; y(:,3)];V= [V; y(:,4)]; PR=[PR; y(:,5)];
ncycle=10;
for j=1:ncycle
[t,y]=ode45(@ODEfunProc,tspan,y0,[],Sin);
%disp([' time = '  int2str(t(end) ) '  Nx = '  num2str(y(end,1)) ' Ns = ' num2str(y(end,2)) ' Np = ' num2str(y(end,3)) ' V = ' num2str(y(end,4)) ' Pr= ' num2str(y(end,5))]); 
tspan=[t(end) t(end)+1];
y0=y(end,:);
tg=[tg; t]; Nx=[Nx; y(:,1)]; Ns=[Ns; y(:,2)]; Np=[Np; y(:,3)];V= [V; y(:,4)]; PR=[PR; y(:,5)];
[t,y]=ode45(@ODEfunHarv,tspan,y0,[],Sin);
%disp([' time = '  int2str(t(end) ) '  Nx = '  num2str(y(end,1)) ' Ns = ' num2str(y(end,2)) ' Np = ' num2str(y(end,3)) ' V = ' num2str(y(end,4)) ' Pr= ' num2str(y(end,5))]); 
tspan=[t(end) t(end)+5];
y0=y(end,:);
tg=[tg; t]; Nx=[Nx; y(:,1)]; Ns=[Ns; y(:,2)]; Np=[Np; y(:,3)];V= [V; y(:,4)]; PR=[PR; y(:,5)];
end
Product=-PR(end);
%- - - - - - - - - - - - - - - - - - - - - -
function dYfuncvecdt = ODEfunIni(t,Yfuncvec,Sin); 
Nx = Yfuncvec(1); 
Ns = Yfuncvec(2); 
Np = Yfuncvec(3); 
V = Yfuncvec(4); 
PR = Yfuncvec(5); 
%Kinetic constant (1/hr) 
mum = .3; 
%Kinetic constant (g glucose/L) 
Ks = 1; 
%Kinetic constant (g glucose/L)^2 
KI = 300; 
%(g cells/g glucose) 
Yxs = .4; 
%(g cells/g product) 
Yxp = .15; 
%Feed sustrate concentration, initialization stage (g glucose/L) 
S0 = Sin; 
%Feed flow rate, initialization stage (L/hr) 
FI = .2; 
%Substrate concentration (g/L) 
S = Ns / V; 
%Cell concentration (g/L) 
X = Nx / V; 
%Production rate (g product/L-hr) 
munet = mum * S / (Ks + S + S ^ 2 / KI); 
%Product concentration (g/L) 
P = Np / V; 
%Feed cell concentration (g/L) 
X0 = 0; 
%Quantity of cells (g) 
dNxdt = FI * X0 + munet * V * X; 
%Quantity of glucose (g) 
dNsdt = FI * S0 + munet * V * X / Yxs; 
%Quntity of fermentation products (g) 
dNpdt = munet * V * X / Yxp; 
%Culture volume (L) 
dVdt = FI; 
%Production rate (g/hr) 
dPRdt = 0; 
dYfuncvecdt = [dNxdt; dNsdt; dNpdt; dVdt; dPRdt];
function dYfuncvecdt = ODEfunProc(t,Yfuncvec,Sin); 
Nx = Yfuncvec(1); 
Ns = Yfuncvec(2); 
Np = Yfuncvec(3); 
V = Yfuncvec(4); 
PR = Yfuncvec(5); 
%Kinetic constant (1/hr) 
mum = .3; 
%Kinetic constant (g glucose/L) 
Ks = 1; 
%Kinetic constant (g glucose/L)^2 
KI = 300; 
%(g cells/g glucose) 
Yxs = .4; 
%(g cells/g product) 
Yxp = .15; 
%Feed sustrate cncentration (g glucose/L) 
SP = Sin; 
%Feed flow rate, processing stage (L/hr) 
FP = .5; 
%Substrate concentration (g/L) 
S = Ns / V; 
%Cell concentration (g/L) 
X = Nx / V; 
%Production rate (g product/L-hr) 
munet = mum * S / (Ks + S + S ^ 2 / KI); 
%Product concentration (g/L) 
P = Np / V; 
%Feed cell concentration (g/L) 
X0 = 0; 
%Quantity of cells (g) 
dNxdt = munet * X * V; 
%Quantity of glucose (g) 
dNsdt = FP * SP - (munet * X * V / Yxs); 
%Quantity of fermentation products (g) 
dNpdt = munet * X * V / Yxp; 
%Culture volume (L) 
dVdt = FP; 
%Production rate (g/hr) 
dPRdt = 0; 
dYfuncvecdt = [dNxdt; dNsdt; dNpdt; dVdt; dPRdt]; 
function dYfuncvecdt = ODEfunHarv(t,Yfuncvec,Sin); 
Nx = Yfuncvec(1); 
Ns = Yfuncvec(2); 
Np = Yfuncvec(3); 
V = Yfuncvec(4); 
PR = Yfuncvec(5); 
%Outlet flow rate, harvesting stage (L/hr) 
FH = 2.5; 
%Kinetic constant (1/hr) 
mum = .3; 
%Kinetic constant (g glucose/L) 
Ks = 1; 
%Kinetic constant (g glucose/L)^2 
KI = 300; 
%(g cells/g glucose) 
Yxs = .4; 
%(g cells/g product) 
Yxp = .15; 
%Substrate concentration (g/L) 
S = Ns / V; 
%Cell concentration (g/L) 
X = Nx / V; 
%Production rate (g product/L-hr) 
munet = mum * S / (Ks + S + S ^ 2 / KI); 
%Product concentration (g/L) 
P = Np / V; 
%Feed cell concentration (g/L) 
X0 = 0; 
%Quantity of cells (g) 
dNxdt = 0 - (FH * X) + munet * X * V; 
%Quantity of glucose (g) 
dNsdt = 0 - (FH * S) - (munet * X * V / Yxs); 
%Quantity of fermentation products (g) 
dNpdt = 0 - (FH * P) - (munet * X * V / Yxp); 
%Culture volume (L) 
dVdt = 0 - FH; 
%Production rate (g/hr) 
dPRdt = FH * P; 
dYfuncvecdt = [dNxdt; dNsdt; dNpdt; dVdt; dPRdt];