MATLAB Examples

indpensim_run.m

Runs IndPensim

Contents

Copyright

Stephen Goldrick September 2014 University of Manchester, Newcastle University and Perceptive Engineering

Please reference "The Development of an Industrial Scale Fed-Batch Fermentation Simulation", Stepen Goldrick, Andrei Stefen, David Lovett, Gary Montague, Barry Lennox, Submitted to the Journal of Biotechnology % September 2014

All rights reserved. Copyright (c) Newcastle University, University of Manchester and Perceptive Engineering.

close all
clc
clear all

Set up simulation flags

Ctrl_flags.SBC= 1;      % Sequential Batch Control flag (SBC) = 1 -  Operator Control (historical Batch data)
                        % Sequential Batch Control flag (SBC) = 0 -  Sequential Batch Control (Fixed trajectories)
Ctrl_flags.IC = 1;      % Initial Conditions flag (IC) = 1 -  Estimated from historical batch data
                        % Initial Conditions flag (IC) = 0 -  Randomly calculated initial conditions must control using SBC
Ctrl_flags.Inhib = 2;   % Inhibition  flag (Inhib)     = 0 - No inhibition
                        % Inhibition  flag (Inhib)     = 1 - Inhibition DO2, T, pH
                        % Inhibition  flag (Inhib)     = 2 - Inhibition of DO2,T,pH,CO_2_L,PAA and N
Ctrl_flags.Dis = 0;     % Disturbance flag  (Dis)      = 0 - No process disturbances
                        % Disturbance flag  (Dis)      = 1 - In batch fluctuations on mu_P, mu_x, c_s, c_oil,abc,PAA_c,Tcin and O2_in

Ctrl_flags.Faults= 0;   % Fault flag  (Faults)         = 0 - No Faults
                        % Fault flag  (Faults)         = 1 - Aeration rate fault
                        % Fault flag  (Faults)         = 2 - Vessel back pressure  fault
                        % Fault flag  (Faults)         = 3 - Substrate feed rate fault
                        % Fault flag  (Faults)         = 4 - Base flowrate fault
                        % Fault flag  (Faults)         = 5 - Coolant flowrate fault
                        % Fault flag  (Faults)         = 6 - All of the above faults
                        % Fault flag  (Faults)         = 7 - Temperature sensor error
                        % Fault flag  (Faults)         = 8 - pH sensor error
Ctrl_flags.Vis= 1;      % Viscosity flag  - Allows to use the recorded viscosity as a process input
                        % Viscosity flag (Vis)          = 0 - uses simulated viscosity
                        % Viscosity flag (Vis)          = 1 - uses recorded viscosity (from batch records) as input

% Off-line measurement sampling rate and analysis delay
Ctrl_flags.Off_line_m =  24;     % Off-line measurement sampling rate (hours)
Ctrl_flags.Off_line_delay =  4;  % Off-line measurement analysis time delay (hours)
% Plot graphs
Ctrl_flags.plots =  1;           % 0 - No plots 1 - plots
% Select industrial batch to import
Batch_no = 3; % Batch data flag - Select historical batch data as input for simulation (1 - Batch 1, 2-Batch 2...10-Batch 10)

% Import and interpolate industrial batch data
if(Ctrl_flags.IC==1)
    disp('Loading industrial batch data...');
    % Load industrial batch data including key batch characteristics
    [batchFile,color_number,V_sf,PAA_c,N_conc_paa,alpha_kla] = Ind_Batch_data(Batch_no);
    % Determines batch length
    T = getIndBatchLength(batchFile);
    % get the length of the batch. Uses the sugar feedrate channel for reference
    h = 0.2; % sampling rate is 12 mins
    % Reload and interpolate batch data
    Xinterp = loadIndBatch1(batchFile, 0:h:T);
    % Smooths volume profile required for  correct F_{discharge} estimation
    Xinterp.Wt.y =  smoothn( Xinterp.Wt.y,V_sf);
    %  caclulates the F_{discharge} rate
    Wt_d = diff(Xinterp.Wt.y);
    %  caclulates the F_{discharge} rate
    for m =1:1:length(Wt_d)
        if(Wt_d(m)/h < -500) % if there is a large negative spike in volume derivative, then there must have been operator intervention
            Xinterp.F_discharge_cal.y(m,1) = Wt_d(m)*5;
        else
            Xinterp.F_discharge_cal.y(m,1) = 0;
        end
    end
Loading industrial batch data...

Calculates the initial growth rates for batch over the first 70 hours

    [mu_X_CER, mu_p] = full_growth_rates2(Xinterp);
    mu_X_max_CER = max(mu_X_CER(1:4));
    mu_P_max = max(mu_p(100:1:350));
    intial_conds = 0.5; % Assumes X_0 is constant (g h^-1)
    x0.mux = mu_X_max_CER;% maximum specific growth rate of biomass [h^-1] (calculated from CER)
    x0.mup = mu_P_max; % maximum specific growth rate of penicillin [h^-1]
    Batch_time = 0:h:T;

Initialising simulation

    x0.S = 1;      % substrate concentration [g L^{-1}]
    x0.DO2 = Xinterp.DO2.y(1); % initial dissolved oxygen concentration [mg L^{-1}]
    x0.X = intial_conds; % initial biomass concentration [g L^{-1}]
    x0.P = 0;       % penicillin concentration [g L^{-1}]
    x0.V =  Xinterp.Wt.y(1)/1.1;   % initial  volume [L]
    x0.Wt = Xinterp.Wt.y(1); % Initial  Weight [Kg]
    x0.CO2 = 0.038;   % initial carbon dioxide concentration offgas  [% of offgas]
    x0.O2 = Xinterp.O2.y(10)/100;% initial oxygen concentration offgas [% of offgas]
    x0.pH = Xinterp.pH.y(1);    % initial pH
    x0.T =Xinterp.T.y(1) ;     % temperature [K]
    x0.a0 = intial_conds*(1/3);   % type a0 biomass concentration [g L^{-1}]
    x0.a1 = intial_conds*(2/3);   % type a1 biomass concentration [g L^{-1}]
    x0.a3 =  0    ;   % type a3 biomass concentration [g L^{-1}]
    x0.a4 = 0 ;   % type a4 biomass concentration [g L^{-1}]
    x0.Culture_age =0; % initial culture age time A_t [hr]
    x0.PAA = Xinterp.PAA.y(1);% initial PAA conc [mg L^{-1}]
    x0.NH3 = Xinterp.NH3.y(1); % initial Nitrogen concentration [mg L^{-1}]

Temperature and pH Set points for Batch

    Ctrl_flags.T_sp =  mean(Xinterp.T.y); % Temperature Set-point (K)
    Ctrl_flags.pH_sp = mean(Xinterp.pH.y); % pH Set-point (-)
end

Standard batch simulation with randomised initial conditions and batch

if(Ctrl_flags.IC==0)
    Ctrl_flags.SBC= 0; % Standard batch simulation must be controlled through sequential batch control
    Ctrl_flags.Vis= 0; % Standard batch simulation must use simulated viscosity
    color_number =1;
    Optimum_Batch_lenght = 250; % Standard batch length (h)
    Batch_length_variation = 25*randn(1);  % Batch length deviation (h)
    T = Optimum_Batch_lenght+Batch_length_variation; % Batch simulation time  (h)
    T= round(T);
    intial_conds = 0.5+0.05*randn; % X_0 (g h^-1)
    x0.mux = 0.41+0.025*randn;% Maximum specific growth rate of biomass [h^-1]
    x0.mup = 0.041+ 0.0025*randn;  % maximum specific growth rate of penicillin [h^-1]
    h = 0.2; % Simulation sampling rate is 12 mins (h)

Initialising simulation

    x0.S = 1 + 0.1*randn;      % substrate concentration [g L^{-1}]
    x0.DO2 = 15 + 0.5*randn  ; % initial dissolved oxygen concentration [mg L^{-1}]
    x0.X = intial_conds + 0.1*randn; % initial biomass concentration [g L^{-1}]
    x0.P = 0;       % penicillin concentration [g L^{-1}]
    x0.V =  5.800e+04 + 500*randn;   % initial  volume [L]
    x0.Wt = 6.2e+04 + 500*randn; % Initial  Weight [Kg]
    x0.CO2 = 0.038 + 0.001*randn;   % initial carbon dioxide concentration offgas  [% of offgas]
    x0.O2 = 0.20+ 0.05*randn;% initial oxygen concentration offgas [% of offgas]
    x0.pH = 6.5 + 0.1*randn;    % initial pH
    x0.T =297+ 0.5*randn;     % temperature [K]
    x0.a0 = intial_conds*(1/3);   % type a0 biomass concentration [g L^{-1}]
    x0.a1 = intial_conds*(2/3);   % type a1 biomass concentration [g L^{-1}]
    x0.a3 =  0    ;   % type a3 biomass concentration [g L^{-1}]
    x0.a4 = 0 ;   % type a4 biomass concentration [g L^{-1}]
    x0.Culture_age =0; % initial culture age time A_t [hr]
    x0.PAA = 1400+50*randn;% initial PAA conc [mg L^{-1}]
    x0.NH3 = 1700+50*randn; % initial Nitrogen concentration [mg L^{-1}]
    alpha_kla = 85 + 10*randn; % kla_constant [-]
    PAA_c =530000  + 20000*randn ; % Concentration of PAA in F_{PAA} [mg/L]
    N_conc_paa = 75000 + 2000*randn;%  Concentration of N in F_{PAA} [mg/L]
    Batch_time = 0:h:T;

Temperature and pH Set points for Batch

    Ctrl_flags.T_sp = 298; % Temperature Set-point (K)
    Ctrl_flags.pH_sp = 6.5; % pH Set-point (-)
end

Creates process disturbances on growth rates as well as process inputs

using a low pass filter

b1 = 1 - 0.995;
a1 = [1 -0.995];
% Penicillin specific growth rate disturbance: with SD of +/- 0.0009 [hr^{-1}]
v = randn(T/h+1, 1);
distMuP = filter(b1,a1,0.03*v);
Xinterp.distMuP = createChannel('Penicillin specific growth rate disturbance','g/Lh','h',Batch_time,distMuP);
% Biomass specific growth rate disturbance: with SD  +/- 0.011 [hr^{-1}]
v = randn(T/h+1, 1);
distMuX = filter(b1,a1,0.25*v);
Xinterp.distMuX = createChannel('Biomass specific  growth rate disturbance','hr^{-1}','h',Batch_time,distMuX);
% Substrate inlet concentration disturbance: +/- 15 [g L^{-1}]
v = randn(T/h+1, 1);
distcs = filter(b1,a1,300*v);
Xinterp.distcs = createChannel('Substrate concentration disturbance ',' g L^{-1}','h',Batch_time,distcs);
% Oil inlet concentration disturbance: +/- 15 [g L^{-1}]
v = randn(T/h+1, 1);
distcoil = filter(b1,a1,300*v);
Xinterp.distcoil = createChannel('Substrate concentration disturbance ',' g L^{-1}','h',Batch_time,distcoil);
% Acid/Base molar inlet concentration disturbance: +/- 0.004 [mol L^{-1}]
v = randn(T/h+1, 1);
distabc = filter(b1,a1,0.2*v);
Xinterp.distabc = createChannel('Acid/Base concentration disturbance ',' g L^{-1}','h',Batch_time,distabc);
v = randn(T/h+1, 1);
distPAA = filter(b1,a1,300000*v);
Xinterp.distPAA = createChannel('Phenylacetic acid concentration disturbance ',' g L^{-1}','h',Batch_time,distPAA);
% PAA inlet concentration disturbance: +/- 120  [g  L^{-1}]
v = randn(T/h+1, 1);
distPAA = filter(b1,a1,300000*v);
Xinterp.distPAA = createChannel('Phenylacetic acid concentration disturbance ',' g L^{-1}','h',Batch_time,distPAA);
% Coolant temperature inlet concentration disturbance: +/- 3 [K]
v = randn(T/h+1, 1);
distTcin = filter(b1,a1,100*v);
Xinterp.distTcin = createChannel('Coolant inlet temperature disturbance ','K','h',Batch_time,distTcin);
% Oxygen inlet concentration: +/- 0.009 [%]
v = randn(T/h+1, 1);
distO_2in = filter(b1,a1,0.02*v);
Xinterp.distO_2in = createChannel('Oxygen inlet concentration','%','h',Batch_time,distO_2in);

Import parameter list

par = Parameter_list(x0,alpha_kla,N_conc_paa,PAA_c);

Runs simulation

disp('Running IndPenSim...');
[Xref] = indpensim(@fctrl_indpensim, Xinterp, x0, h, T,2,par,Ctrl_flags);
Running IndPenSim...

Plots graphs

disp('Plotting graphs...');
if Ctrl_flags.plots  ==1
    if   Ctrl_flags.SBC == 0
        cmap = lines(11);
        figure(1)
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        set(gcf,'name','Biomass')
        hold on
        plot(Xref.X.t,Xref.X.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Biomass concentrations-{\itX},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('X sim')





        figure(100)
        set(gcf,'name','Biomass regions')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.a0.t, Xref.a0.y,'r')
        plot(Xref.a1.t, Xref.a1.y,'y')
        plot(Xref.a3.t, Xref.a3.y,'k')
        plot(Xref.a4.t, Xref.a4.y,'b')
        plot(Xref.X.t,Xref.X.y,'m')
        total_X = Xref.a0.y +  Xref.a1.y+ Xref.a3.y +Xref.a4.y;
        plot(Xref.a4.t,total_X,'g')
        set(get(gca,'YLabel'),'String','Biomass concentrations-{\itX,A_0,A_1,A_3,A_4},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('a0','a1','a3','a4','X','Total X','Location','Best');



        figure(2)
        set(gcf,'name','Penicillin')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.P.t,Xref.P.y,'Color',cmap(color_number,:))
        plot(Xref.P_offline.t,Xref.P_offline.y,'v-.b','linewidth',1.5,'MarkerEdgeColor','b');
        set(get(gca,'YLabel'),'String','Penicillin concentration-{\itP},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('P sim','P  sim off-line','Location','Best')

        figure(3)
        hold on
        set(gcf,'name','Substrate feed rates')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        [AX,H1,H2] = plotyy(Xref.Fs.t,Xref.Fs.y,Xref.Foil.t,Xref.Foil.y,'plot');
        set(get(AX(1),'Ylabel'),'String','Slow Decay')
        set(H1,'Color',cmap(color_number,:))
        set(H2,'Color',cmap(color_number+1,:))
        xlabel('Batch time (hr)')
        legend('Sugar flow rate F_s', 'Soybean flow rate F_{oil}','Location','Best')
        set(get(gca,'YLabel'),'String','Sugar and Soybean flow rate-{\itF_{s},F_{oil}},(L h^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')




        figure(40)
        hold on
        set(gcf,'name','Substrate conc')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        plot(Xref.S.t,Xref.S.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Substrate concentration-{\its},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('substrate sim','Location','Best')


        figure(4)
        hold on
        set(gcf,'name','DO_2')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        plot(Xref.DO2.t, Xref.DO2.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Dissolved oxygen concentration-{\itDO_2},(mg L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        % Estimated the DO_2 limit in vesssel
        % Calculating liquid height in vessel
        h_b = (Xref.V.y./1000)/(pi()*(2.1^2));
        h_b = h_b*(1-0.1); % (ungassed height)
        % Calculating log mean pressure of vessel [bar]
        pressure = Xref.pressure.y;
        pho_b = 1100;
        pressure_bottom  =  1+ pressure + ((pho_b*h_b)*9.81*10^(-5)); % [bar]
        pressure_top = 1+ pressure; % [bar]
        log_mean_pressure = (pressure_bottom - pressure_top)./(log(pressure_bottom./pressure_top));
        total_pressure = log_mean_pressure;
        hold on
        Henrys_c = 0.021;
        DOstar_tp_3 =  0.3*(((total_pressure)*0.21)/Henrys_c);%in mg/l (Henry's constant has the units of (bar L mg-1)
        plot(Xref.CO2_d.t,DOstar_tp_3)
        legend('DO_2 sim','DO_2 critical limit (30%)','Location','Best')



        figure(5)
        hold on
        set(gcf,'name','Fg and ')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        [AX,H1,H2] = plotyy(Xref.Fg.t, Xref.Fg.y,Xref.Fg.t, Xref.pressure.y,'plot');
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        set(get(AX(1),'Ylabel'),'String','Aeration rate, -{\itF_g},(m^3 min^{-1})')
        set(get(AX(2),'Ylabel'),'String','vessel back pressure - {\itP_1} (bar)')
        legend('Air Flow', 'Pressure','Location','Best')








        figure(6)
        set(gcf,'name','Temp')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.T.t, Xref.T.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Temperature- {\it T_b},(K)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('T sim','Location','Best')

        figure(7)
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        set(gcf,'name','Cooling water flow')
        set(get(gca,'YLabel'),'String','Cooling water-{\itF_c},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        hold on
        plot(Xref.Fc.t,Xref.Fc.y,'Color',cmap(color_number,:))
        legend('Cooling water flow','Location','Best')

        figure(8)
        set(gcf,'name','hot water  flow')
         set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fh.t,Xref.Fh.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Heating water-{\itF_h},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Heating water flow','Location','Best')

        figure(9)
        set(gcf,'name','pH')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.pH.t, Xref.pH.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','pH,')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('pH sim','Location','Best')

        figure(10)
        set(gcf,'name','Fa')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fa.t,Xref.Fa.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Acid flow rate -{\itF_a},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Acid  flow','Location','Best')
        figure(11)
        set(gcf,'name','Fb')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fb.t,Xref.Fb.y,'Color',cmap(color_number,:))
         set(get(gca,'YLabel'),'String','Base flow rate -{\itF_b},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Base flow rate','Location','Best')

        figure(12)
        set(gcf,'name','Viscosity')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Viscosity.t, Xref.Viscosity.y,'Color',cmap(color_number,:))
        plot(Xref.Viscosity_offline.t,Xref.Viscosity_offline.y,'v-.b','linewidth',1.5,'MarkerEdgeColor','b');
        set(get(gca,'YLabel'),'String','Viscosity - {\itu_{app}},(cP)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('viscosity sim',' viscosity off-line sim','Location','Best')

        figure(13)
        set(gcf,'name','Water for injection')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fw.t, Xref.Fw.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Water for injection - {\itF_{w}},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('water for injection','Location','Best')

        figure(14)
        set(gcf,'name','PAA')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.PAA.t, Xref.PAA.y,'Color',cmap(color_number,:))
        plot(Xref.PAA_offline.t,Xref.PAA_offline.y,'v-.b','linewidth',1.5,'MarkerEdgeColor','b');
        set(get(gca,'YLabel'),'String','Phenylacetic acid conc - {\itPAA},(mg L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('PAA', 'PAA off-line','Location','Best')




        figure(15)
        set(gcf,'name','Fpaa')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fpaa.t, Xref.Fpaa.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Phenylacetic acid flow rate - {\itF_{PAA}},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Fpaa','Location','Best')

        figure(16)
        set(gcf,'name','NH3')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.NH3.t, Xref.NH3.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Nitrogen conc - {\itN},(mg L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('NH_3 sim','Location','Best')

        figure(17)
        set(gcf,'name','CO2')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.CO2outgas.t, (Xref.CO2outgas.y),'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','CO_2 out-gas - {\itCO_{2,out}},(%)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
         legend('CO_{2,out} sim','Location','Best')



        figure(18)
        set(gcf,'name','O2 outgas')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.O2.t, (Xref.O2.y*100),'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','O_2 out-gas - {\itO_{2,out}},(%)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
         legend('O_{2,out} sim','Location','Best')


        figure(19)
        set(gcf,'name','Volume')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.V.t, Xref.V.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Volume - {\itV},(L)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Volume sim','Location','Best')

        figure(20)
        set(gcf,'name','Weight')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Wt.t, Xref.Wt.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Weight - {\itWT},(kg)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Weight sim','Location','Best')

         figure(24)
        set(gcf,'name','OUR')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.OUR.t,abs(Xref.OUR.y),'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Oxygen uptake rate  - {\itOUR},(g/min)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('OUR sim','Location','Best')


        figure(25)
        set(gcf,'name','CO_2 dissolved')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.CO2_d.t,Xref.CO2_d.y*1000,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Dissolved carbon dioxide - {\itCO_{2,L}},(-)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('CO_2 dissolved sim','Location','Best')

        figure(26)
        set(gcf,'name','CER')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.CER.t,Xref.CER.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Carbon evolution rate - {\itCER},(g/min)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('CER sim')

    end


    if (Ctrl_flags.SBC==1)
        % load Batch_5
        cmap = lines(12);
        % plotting the main results:
        figure(1)
        set(gcf,'name','Biomass')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.X_CER.t,Xinterp.X_CER.y,'Color',cmap(color_number+1,:))
        plot(Xref.X.t,Xref.X.y,'Color',cmap(color_number,:))
        total_Biomass  = Xref.a0.y+Xref.a1.y+Xref.a3.y+Xref.a4.y;
        %plot(Xref.DO2.t,total_Biomass)
        set(get(gca,'YLabel'),'String','Biomass concentrations-{\itX},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('X softsensor', 'X sim')
        figure(100)
        set(gcf,'name','Biomass regions')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.a0.t, Xref.a0.y,'r')
        plot(Xref.a1.t, Xref.a1.y,'y')
        plot(Xref.a3.t, Xref.a3.y,'k')
        plot(Xref.a4.t, Xref.a4.y,'b')
        plot(Xref.X.t,Xref.X.y,'m')
        total_X = Xref.a0.y +  Xref.a1.y+ Xref.a3.y +Xref.a4.y;
        plot(Xref.a4.t,total_X,'g')
        set(get(gca,'YLabel'),'String','Biomass concentrations-{\itX,A_0,A_1,A_3,A_4},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('a0','a1','a3','a4','X','Total X');


        figure(2)
        set(gcf,'name','Penicillin')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.P.t,Xinterp.P.y,'Color',cmap(color_number+1,:))
        plot(Xinterp.P_offline.t,Xinterp.P_offline.y,'v-.r','linewidth',1.5,'MarkerEdgeColor','r');
        plot(Xref.P.t,Xref.P.y,'Color',cmap(color_number,:))
        plot(Xref.P_offline.t,Xref.P_offline.y,'v-.b','linewidth',1.5,'MarkerEdgeColor','b');
        set(get(gca,'YLabel'),'String','Penicillin concentration-{\itP},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('P ','P off-line','P sim ','P sim off-line','Location','Best')

        figure(3)
        hold on
        set(gcf,'name','Substrate flow rate')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        plot(Xref.Fs.t,Xref.Fs.y,'Color',cmap(color_number,:))
        plot(Xref.Foil.t,Xref.Foil.y,'Color',cmap(color_number+2,:))
        set(get(gca,'YLabel'),'String','Sugar and Soybean flow rate-{\itF_{s},F_{oil}},(L h^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Sugar flow rate F_s', 'Soybean flow rate F_{oil}','Location','Best')


        figure(40)
        hold on
        set(gcf,'name','Substrate conc')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        plot(Xref.S.t,Xref.S.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Substrate concentration-{\its},(g L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('substrate sim','Location','Best')
        legend('Substrate conc')


        figure(4)
        hold on
        set(gcf,'name','DO_2')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        plot(Xinterp.DO2.t, Xinterp.DO2.y,'Color',cmap(color_number+1,:))
        plot(Xref.DO2.t, Xref.DO2.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Dissolved oxygen concentration-{\itDO_2},(mg L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        % Estimated the DO_2 limit in vessel
        % Calculating liquid height in vessel
        h_b = (Xref.V.y./1000)/(pi()*(2.1^2));
        h_b = h_b*(1-0.1); % (ungassed height)
        % Calculating log mean pressure of vessel [bar]
        pressure = Xref.pressure.y;
        pho_b = 1100;
        pressure_bottom  =  1+ pressure + ((pho_b*h_b)*9.81*10^(-5)); % [bar]
        pressure_top = 1+ pressure; % [bar]
        log_mean_pressure = (pressure_bottom - pressure_top)./(log(pressure_bottom./pressure_top));
        total_pressure = log_mean_pressure;
        hold on
        Henrys_c = 0.021;
        DOstar_tp_3 =  0.3*(((total_pressure)*0.21)/Henrys_c);%in mg/l (henrys constant has the units of (bar L mg-1)
        plot(Xref.CO2_d.t,DOstar_tp_3)
        legend('DO_2', 'DO_2 sim','DO_2 critical limit (30%)','Location','Best')




        figure(5)
        hold on
        set(gcf,'name','Fg and ')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        [AX,H1,H2] = plotyy(Xref.Fg.t, Xref.Fg.y,Xref.Fg.t, Xref.pressure.y,'plot');
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        set(get(AX(1),'Ylabel'),'String','Aeration rate, -{\itF_g},(m^3 min^{-1})')
        set(get(AX(2),'Ylabel'),'String','vessel back pressure - {\itP_1} (bar)')
        legend('Air Flow', 'Pressure','Location','Best')



        figure(6)
        set(gcf,'name','Temp')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.T.t, Xinterp.T.y,'Color',cmap(color_number+1,:))
        plot(Xref.T.t, Xref.T.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Temperature- {\it T_b},(K)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('T', 'T sim','Location','Best')

        figure(7)
        set(gcf,'name','Cooling water flow')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fc.t,Xref.Fc.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Cooling water-{\itF_c},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
       legend('Cooling water flow','Location','Best')


        figure(8)
        set(gcf,'name','hot water  flow')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fh.t,Xref.Fh.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Heating water-{\itF_h},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Heating water flow','Location','Best')


        figure(9)
        set(gcf,'name','pH')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.pH.t, Xinterp.pH.y,'Color',cmap(color_number+1,:))
        plot(Xref.pH.t, Xref.pH.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','pH,')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('pH', 'pH sim','Location','Best')

        figure(10)
        set(gcf,'name','Fa')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fa.t,Xref.Fa.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Acid flow rate -{\itF_a},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Acid  flow','Location','Best')

        figure(11)
        set(gcf,'name','Fb')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fb.t,Xref.Fb.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Base flow rate -{\itF_b},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Base flow rate','Location','Best')



        figure(12)
        set(gcf,'name','Viscosity')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.viscosity.t, Xinterp.viscosity.y,'Color',cmap(color_number+1,:))
        plot(Xinterp.viscosity_offline.t,Xinterp.viscosity_offline.y,'v-.r','linewidth',1.5,'MarkerEdgeColor','r');
        plot(Xref.Viscosity.t, Xref.Viscosity.y,'Color',cmap(color_number,:))
        plot(Xref.Viscosity_offline.t,Xref.Viscosity_offline.y,'v-.b','linewidth',1.5,'MarkerEdgeColor','b');
        set(get(gca,'YLabel'),'String','Viscosity - {\itu_{app}},(cP)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Viscosity','Viscosity offline','Viscosity sim','Viscosity offline sim','Location','Best')


        figure(13)
        set(gcf,'name','Water for injection')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.Fw.t, Xref.Fw.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Water for injection - {\itF_{w}},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('water for injection','Location','Best')
        legend('water for injection')

        figure(14)
        set(gcf,'name','PAA')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.PAA.t, (Xinterp.PAA.y),'Color',cmap(color_number+1,:))
        plot(Xinterp.PAA_offline.t,Xinterp.PAA_offline.y,'v-.r','linewidth',1.5,'MarkerEdgeColor','r');
        plot(Xref.PAA.t, Xref.PAA.y,'Color',cmap(color_number,:))
        plot(Xref.PAA_offline.t,Xref.PAA_offline.y,'v-.b','linewidth',1.5,'MarkerEdgeColor','b');
        set(get(gca,'YLabel'),'String','Phenylacetic acid conc - {\itPAA},(mg L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('PAA','PAA offline','PAA sim','PAA sim off')

        figure(15)
        set(gcf,'name','Fpaa')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.Fpaa.t, Xinterp.Fpaa.y,'Color',cmap(color_number+1,:))
        plot(Xref.Fpaa.t, Xref.Fpaa.y,'Color',cmap(color_number,:))
         set(get(gca,'YLabel'),'String','Phenylacetic acid flow rate - {\itF_{PAA}},(L hr^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Fpaa','Location','Best')

        figure(16)
        set(gcf,'name','NH3')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.NH3.t, Xinterp.NH3.y,'Color',cmap(color_number+1,:))
        plot(Xinterp.NH3_offline.t,Xinterp.NH3_offline.y,'v-.r','linewidth',1.5,'MarkerEdgeColor','r');
        plot(Xref.NH3.t, Xref.NH3.y,'Color',cmap(color_number,:))
        plot(Xref.NH3_offline.t,Xref.NH3_offline.y,'v-.b','linewidth',1.5,'MarkerEdgeColor','b');
        set(get(gca,'YLabel'),'String','Nitrogen conc - {\itN},(mg L^{-1})')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('NH3', 'NH3 off-line','NH_3 sim','NH_3 sim offline','Location','Best')


        figure(17)
        set(gcf,'name','CO2')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.CO2.t, Xinterp.CO2.y,'Color',cmap(color_number+1,:))
        plot(Xref.CO2outgas.t, (Xref.CO2outgas.y),'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','CO_2 out-gas - {\itCO_{2,out}},(%)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
         legend('CO_{2,out}','CO_{2,out} sim','Location','Best')

        figure(18)
        set(gcf,'name','O2 outgas')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.O2.t, Xinterp.O2.y,'Color',cmap(color_number+1,:))
        plot(Xref.O2.t, (Xref.O2.y*100),'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','O_2 out-gas - {\itO_{2,out}},(%)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('O_{2,out}','O_{2,out} sim','Location','Best')
        legend('OUR', 'OUR simulated')



        figure(19)
        set(gcf,'name','Volume')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.V.t, Xref.V.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Volume - {\itV},(L)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Volume sim','Location','Best')


        figure(20)
        set(gcf,'name','Weight')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.Wt.t, Xinterp.Wt.y,'Color',cmap(color_number+1,:))
        plot(Xref.Wt.t, Xref.Wt.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Weight - {\itWT},(kg)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('Weight', 'Weight simulated','Location','Best')



        figure(24)
        set(gcf,'name','OUR')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.OUR.t, Xinterp.OUR.y,'Color',cmap(color_number+1,:))
        plot(Xref.OUR.t,abs(Xref.OUR.y),'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Oxygen uptake rate  - {\itOUR},(g/min)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('OUR','OUR sim','Location','Best')


        figure(25)
        set(gcf,'name','CO_2 dissolved')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xref.CO2_d.t,Xref.CO2_d.y*1000,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Dissolved Carbon dioxide- {\itCO_{2,L}},(mg/L)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('CO_2 dissolved')


         figure(26)
        set(gcf,'name','CER')
        set(gcf, 'color', [1 1 1])
        set(gca,'color', [1 1 1])
        hold on
        plot(Xinterp.CER.t, Xinterp.CER.y,'Color',cmap(color_number+1,:))
        plot(Xref.CER.t,Xref.CER.y,'Color',cmap(color_number,:))
        set(get(gca,'YLabel'),'String','Carbon evolution rate - {\itCER},(g/min)')
        set(get(gca,'XLabel'),'String',' Batch time (hr) ')
        legend('CER','CER sim')

    end
end
Plotting graphs...