bioreactor_consOneD.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DFBAlab: Dynamic Flux Balance Analysis laboratory                       %
% Process Systems Engineering Laboratory, Cambridge, MA, USA              %
% July 2014                                                               %
% Written by Jose A. Gomez and Kai Höffner                                %
%                                                                         % 
% This code can only be used for academic purposes. When using this code  %
% please cite:                                                            %
%                                                                         %
% Gomez, J.A., Höffner, K. and Barton, P. I.                              %
% DFBAlab: A fast and reliable MATLAB code for Dynamic Flux Balance       %
% Analysis. Submitted.                                                    %
%                                                                         %
% COPYRIGHT (C) 2014 MASSACHUSETTS INSTITUTE OF TECHNOLOGY                %
%                                                                         %
% Read the LICENSE.txt file for more details.                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Gas fermentation in a bubble column reactor - main program
% Written by Jin Chen and Michael A. Henson
% September 2017
function new_output = bioreactor_consOneD(tr,ug)   %[Temp,ug]

addpath(genpath('/home/kudva.7/Desktop/PaulsonLab/Akshay/utilities'))

% Changed the input sequence for the maximin case
%tr = input(1,1);
%klac = input(1,2);
%ug = input(1,2);
%ul = input(1,4);

plot_figure = 0;
% Set number of models and simulation time
%tic    
nmodel = 50;
tspan = [0,1000];
N = nmodel+1;                   % number of reactor discretization points
ns = 4;                         % number of state variables change with location

% Load models. These should be .mat files generated by the COBRA toolbox. 
% When generating these files using the COBRA toolbox, a big number is used
% as infinity. This number should be fed to the DB vector (Default bound).
% INFO.DB = DB;

for i = 1:nmodel
    load modelwt.mat           % this GEM is from www.sciencedirect.com/science/article/pii/S2405471217301400
    model{i} = modelwt; 
    DB(i) = 1000;
end

%% exID array
% You can either search the reaction names by name or provide them directly
% in the exID array.
% RxnNames = {'EX_glc(e)', 'EX_ac(e)', 'biomass'};
% for i = 1:length(RxnNames)
%    [a,exID(i)] = ismember(RxnNames(i),model.rxns);
% end
% INFO.exID = exID;
% Lower bounds and upper bounds for these reactions should be provided in
% the RHS code. 

for i = 1:nmodel
    exID{i}=[1026, 960, 854, 972, 969, 986, 658, 51];  
end

% This codes solves the LPs in standard form. Bounds on exchange fluxes in 
% the exID array can be modified directly on the first 2*n rows where n is 
% the number of exchange fluxes. Order will be lower bound, upper bound, 
% lower bound, upper bound in the same order as exID. 
%
% NOTE: All bounds on fluxes in the exID arrays are relaxed to -Inf and 
% + Inf. These bounds need to be updated if needed in the RHS file.

%% Cost vectors
% Usually the first cost vector will be biomass maximization, but it can
% be any other objective. The CPLEX objects will minimize by default. 
% Report only nonzero elements. 
% The structure should be:
% C{model} = struct
% Element C{k}(i) of C, is the cost structure i for model k. 
% C{k}(i).sense = +1 for minimize, or -1 for maximize.
% C{k}(i).rxns = array containing the reactions in this objective. 
% C{k}(i).wts = array containing coefficients for reactions reported in 
% rxns. Both arrays should have the same length. 
% Example, if:
% C{k}(i).rxns = [144, 832, 931];
% C{k}(i).wts = [3, 1, -1];
% Then the cost vector for this LP will be:
% Cost{k}(i) = 3*v_144 + v_832 - v_931 (fluxes for model k). 
% This cost vector will be either maximized or minimized depending on the
% value of C{k}(i).sense.

% In SBML files, usually production fluxes are positive and uptake fluxes
% are negative. Keep in mind that maximizing a negative flux implies 
% minimizing its absolute value.
% Different models can have different number of objectives. 
% INFO.C = C; 

minim = 1;
maxim = -1;

for i=1:nmodel
    
    % Maximize growth
    C{i}(1).sense = maxim;
    C{i}(1).rxns = [1026];
    C{i}(1).wts = [1];
    % Maximize CO uptake
    C{i}(2).sense = minim;
    C{i}(2).rxns = [960];
    C{i}(2).wts = [1];
    % Minimize CO2 production
    C{i}(3).sense = minim;
    C{i}(3).rxns = [854];
    C{i}(3).wts = [1];
    % Minimize acetate uptake
    C{i}(4).sense = minim;
    C{i}(4).rxns = [972];
    C{i}(4).wts = [1];
    % Minimize ethanol uptake
    C{i}(5).sense = minim;
    C{i}(5).rxns = [969];
    C{i}(5).wts = [1];
    % Minimize 2,3-Butanediol uptake
    C{i}(6).sense = minim;
    C{i}(6).rxns = [986];
    C{i}(6).wts = [1];
    % Minimize lactate uptake
    C{i}(7).sense = minim;
    C{i}(7).rxns = [658];
    C{i}(7).wts = [1];

end

%%
% Pass parameters

INFO.nmodel = nmodel;
INFO.DB = DB;
INFO.exID = exID;
INFO.C = C; 

% Set operating conditions
L = 1.06;                         % reactor length, m 
%ug = 12.3;                        % 12.3 constant gas velocity, m/h
ul = -158.5;                      % -158.5 Liquid recyle rate, m/h
D = 0.06;                         % dillution rate,1/h
eg = 0.309;                       % gas holdup, %
Area = 0.002436;                  % reactor cross-sectional area, m^2
zs = L/(N-1);                     % spatial step size, m
dl = 4.5;                         % liquid dispersion coefficient, 0.25 m^2/h
%tr = 310.15;                      %310.15; % temperature,K
pL = 1.013e5;                     % pressure at top of column, Pa
pc = 0.5;                         % partial pressure of CO %(60+2*(mm-1))/100;
Hc = 8.0e-4;                      % Henry's constant for CO in water, mol/L*atm
pc2 = 0.2;                        % partial pressure of CO2
Hc2 = 2.5e-2;                     % Henry's constant for CO2 in water, mol/L*atm
Qmedia = L*Area*D;

% Set mass transfer coefficients
klac = 523;                       % 523 CO gas-liquid mass transfer coefficient, h-1
klac2 = klac;                     % CO2 gas-liquid mass transfer coefficient, h-1

% Calculate gas and liquid holdups
el = 1-eg;
po = pL+1000*9.81*L*el;
cgi = pc*po/8.314/tr;           % CO feed concentration, mol/m3 = mmol/L
c2gi = pc2*po/8.314/tr;         % CO2 feed concentration, mol/m3 = mmol/L

% Form condition vector
condit = [klac,Hc,klac2,Hc2,tr,zs,N,ns,ug,ul,dl,eg,el,cgi,c2gi,Area,Qmedia,D];

% Set uptake parameters
vcm = 35;  
Kmc = 0.02;
Kic = 0.6;     
param = [
     vcm                 % maximum CO uptake rate, mmol/g/h; 
     Kmc                 % CO uptake saturation constant, mmol/L;
     vcm                 % maximum CO2 uptake rate, mmol/g/h 
     Kmc                 % CO2 uptake saturation constant, mmol/L;                 
     Kic                 % CO inhibition constant, mmol/L
     ];

% Pass parameters

INFO.param = param;
INFO.ns = ns;
INFO.N = N;
INFO.condit = condit;

% Set initial conditions

yo = [];

for i = 1:N
    po = pL + 1000*9.81*zs*(N-i)*el;   
    cgii = pc*po/8.314/tr;           
    c2gii = pc2*po/8.314/tr; 
    clsi = cgii*8.314*tr*Hc*1000/1.013e5;
    c2lsi = c2gii*8.314*tr*Hc2*1000/1.013e5;  
    yz = [cgii c2gii clsi c2lsi];
    yo = [yo yz];
end
yo = [yo, 0.1, 0, 0, 0, 0, 0];
Y0 = yo;

% CPLEX Objects construction parameters

INFO.LPsolver = 1; % CPLEX = 0, Gurobi = 1.
                   % CPLEX works equally fine with both methods.
                   % Gurobi seems to work better with Method = 1, and 
                   % Mosek with Method = 0.
INFO.tol = 1E-9; % Feasibility, optimality and convergence tolerance for Cplex (tol>=1E-9). 
                 % It is recommended it is at least 2 orders of magnitude
                 % tighter than the integrator tolerance. 
                 % If problems with infeasibility messages, tighten this
                 % tolerance.
INFO.tolPh1 = INFO.tol; % Tolerance to determine if a solution to phaseI equals zero.
                   % It is recommended to be the same as INFO.tol. 
INFO.tolevt = 10*INFO.tol; % Tolerance for event detection. Has to be greater 
                   % than INFO.tol.

% You can modify the integration tolerances here.
% If some of the flows become negative after running the simulation once
% you can add the 'Nonnegative' option.

options = odeset('AbsTol',1E-6,'RelTol',1E-6, ...
'Events',@evts,'OutputSel',1,'OutputFcn',@odeplot); % 1 is for tracking the first column

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[model,INFO] = ModelSetupM(model,Y0,INFO);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if INFO.LPsolver == 0
    [INFO] = LexicographicOpt(model,INFO);
elseif INFO.LPsolver == 1
    [INFO] = LexicographicOptG(model,INFO);
else
    display('Solver not currently supported.');
end

tint = 0;
TF = [];
YF = [];
figure(1000)
while tint<tspan(2)
    
% Look at MATLAB documentation if you want to change solver.
% ode15s is more or less accurate for stiff problems. 

    [T,Y] = ode15s(@DRHS,tspan,Y0,options,INFO); 
    TF = [TF;T];
    YF = [YF;Y];
    tint = T(end);
    tspan = [tint,tspan(2)];
    Y0 = Y(end,:);
    if tint == tspan(2)
        break;
    end
    
% Determine model with basis change

    value = evts(tint,Y0,INFO);
    [jjj,j] = min(value);
    ct = 0;
    k = 0;
    while j>ct
        k = k + 1;
        ct = ct + size(model{k}.A,1);
    end
    INFO.flagbasis = k;
    fprintf('Basis change at time %d. \n',tint);
    
% Update b vector 

[INFO] = bupdate(tint,Y0,INFO);

% Perform lexicographic optimization

if INFO.LPsolver == 0
    [INFO] = LexicographicOpt(model,INFO);
elseif INFO.LPsolver == 1
    [INFO] = LexicographicOptG(model,INFO);
else
    display('Solver not currently supported.');
end

end


T = TF;
Y = YF;

coeff=[];
for iii=1:length(T)
    [~, coeff(iii,:)]=DRHS(T(iii),Y(iii,:),INFO);
end
Y(:,(N-1)*ns+3)=coeff(:,1);
Y(:,(N-1)*ns+4)=coeff(:,2);

%% Plotting
if plot_figure
figure(1)
set(gcf,'Position',[0 0 1000 600])
subplot(3,3,1)
hold on
plot(T,Y(:,ns*N+1))
xlabel('Time [h]')
ylabel('Biomass [g/L]')
xlim([0 tspan(2)])

subplot(3,3,2)
hold on
plot(T,Y(:,ns*(N-1)+3))
xlabel('Time [h]')
ylabel('CO in liquid [mmol/L]')
xlim([0 tspan(2)])

subplot(3,3,3)
hold on
plot(T,Y(:,ns*(N-1)+1))
xlabel('Time [h]')
ylabel('CO in gas [mmol/L]')
xlim([0 tspan(2)])

subplot(3,3,4)
hold on
plot(T,Y(:,ns*(N-1)+4))
xlabel('Time [h]')
ylabel('CO2 in liquid [mmol/L]')
xlim([0 tspan(2)])

subplot(3,3,5)
hold on
plot(T,Y(:,ns*(N-1)+2))
xlabel('Time [h]')
ylabel('CO2 in gas [mmol/L]')
xlim([0 tspan(2)])

subplot(3,3,6)
hold on
plot(T,Y(:,ns*N+2))
xlabel('Time [h]')
ylabel('Acetate [g/L]')
xlim([0 tspan(2)])

subplot(3,3,7)
hold on
plot(T,Y(:,ns*N+3))
xlabel('Time [h]')
ylabel('Ethanol [g/L]')
xlim([0 tspan(2)])

subplot(3,3,8)
hold on
plot(T,Y(:,ns*N+4))
xlabel('Time [h]')
ylabel('BDOH [g/L]')
xlim([0 tspan(2)])

subplot(3,3,9)
hold on
plot(T,Y(:,ns*N+5))
xlabel('Time [h]')
ylabel('Lactate [g/L]')
xlim([0 tspan(2)])
end

% Calculate spatial profiles at final time
for i=1:N
    zt(i)=(i-1)*zs;                % zt from 1 to L
    yfin1(i)=Y(end,ns*(i-1)+1);
    yfin2(i)=Y(end,ns*(i-1)+2);
    yfin3(i)=Y(end,ns*(i-1)+3);
    yfin4(i)=Y(end,ns*(i-1)+4);
end
   
%%%%%%%%%%%%%%%%%%%%%%% Edited for Constrained Minmax %%%%%%%%%%%%%%%%%%%%%%%%
% Time to steady state

ethanol_conc = Y(:,ns*N+3);
acetate_conc = Y(:,ns*N+2);
selectivity_traj = ethanol_conc./acetate_conc;

temp_vec = [T,ethanol_conc,acetate_conc,selectivity_traj];

skip_to_index = 400;

size_temp_vec = size(temp_vec,1); % The size of this vector is a variable
temp_vec = [[1:size_temp_vec]', temp_vec];

selectivity_v1 = [temp_vec(skip_to_index:end,1),temp_vec(skip_to_index:end,5)];
selectivity_v2 = [abs(diff(selectivity_v1(:,2)));0];

% index_int = find(selectivity_v2<0.0001);
% p=find(diff(index_int)==1);

output =  temp_vec(end,2:end);


% output(:,3) = [];
% 
% v = output(:, 1);
% output(:, 1) = -output(:, 2);
% output(:, 2) = v;

% FinalResulltVals: ethanol_ss = output(:,1) + 14  ;time_ss = output(:,2) - 810; selectivity = 1.515 -output(:,3) ;
% Relaxation Again
ethanol_ss = output(:,2)  ;time_ss = output(:,1) ; selectivity = output(:,4) ; acetate_ss = output(:,3);

new_output = [ethanol_ss,acetate_ss];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if plot_figure
figure(2)
set(gcf,'Position',[0 0 600 400])
subplot(2,2,1)
hold on
plot(zt,yfin1)
xlabel('Location [m]')
ylabel('CO in gas [mmol/L]')
xlim([0 L])

subplot(2,2,2)
hold on
plot(zt,yfin2)
xlabel('Location [m]')
ylabel('CO2 in gas [mmol/L]')
xlim([0 L])

subplot(2,2,3)
hold on
plot(zt,yfin3)
xlabel('Location [m]')
ylabel('CO in liquid [mmol/L]')
xlim([0 L])

subplot(2,2,4)
hold on
plot(zt,yfin4)
xlabel('Location [m]')
ylabel('CO2 in liquid [mmol/L]')
xlim([0 L])
end
end