import numpy as np
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

#Simulate fed-batch operation
# Specify the simulation time (hrs)
tspan = np.linspace(0,37,371); t=tspan
# Specify values of control variables [Qin0 Qin Xtin Xvin Qe Sin Fc Fair Tin Tcin]
u0 = [0.0,0.0,0.0,0.0,0.0,400,40,60000,30,15]
# Specify initial conditions [Xt Xv S P Oliq Ogas T Tc Vl Gloss MP]
x0 = [0.1,0.1,50,0,0.0065,0.305,30,20,1000,0.0,0.0]
ux0 = (u0 + x0)

#Qin = 15*heaviside(t-5) + 5*heaviside(t-10)
#      - 6*heaviside(t-20) - 14*heaviside(t-35)
Qin = np.zeros_like(tspan)
Qin[np.where(tspan>=5)]  += 15
Qin[np.where(tspan>=10)] += 5
Qin[np.where(tspan>=20)] -= 6
Qin[np.where(tspan>=35)] -= 14
QinInterp = interp1d(tspan,Qin,bounds_error=False)

def ethanol(x,t,Qin0,Qin,Xtin,Xvin,Qe,Sin,Fc,Fair,
            Tin,Tcin,Xt0,Xv0,S0,P0,Oliq0,Ogas0,T0,
            Tc0,Vl0,Sf_cum0,Time0):
    ##    #Initial Conditions
    ##    Xt0 = u[10]     # Initial total cellular biomass, [g L-1]
    ##    Xv0 = u[11]     # Initial viable cellular biomass, [g L-1]
    ##    S0 = u[12]      # Initial substrate/Glucose concentration, [g L-1]
    ##    P0 = u[13]      # Initial product/Ethanol concentration, [g L-1]
    ##    Oliq0 = u[14]   # Initial Dissolved oxygen concentration, [g L-1]
    ##    Ogas0 = u[15]   # Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
    ##    T0 = u[16]      # Initial Temperature in the bioreactor, [oC]
    ##    Tc0 = u[17]     # Initial Temperature of the cooling agent in the jacket, [oC]
    ##    Vl0 = u[18]     # Initial Culture volume in the bioreactor, [L]
    ##    Sf_cum0 = u[19] # Initial Cumulative substrate/glucose fed to the bioreactor, [g]
    ##    Time0 = u[20]   # Initial batch time, [h]
    ##    
    ##    #Control variables
    ##    Qin0 = u[0]        # Volumetric inflow rate, [l/h-1]
    ##    Qin = u[1]         # Volumetric inflow rate, [l/h-1]
    ##    Xtin = u[2]        # Total biomass concentration in the bioreactor feed, [g L-1]
    ##    Xvin = u[3]        # Viable biomass concentration in the bioreactor feed, [g L-1]
    ##    Qe = u[4]          # Volumetric outflow rate, [l/h-1]
    ##    Sin = u[5]         # Substrate/Glucose concentration in bioreactor feed, [g L-1]
    ##    Fc = u[6]          # Cooling agent inlet volumetric flowrate, [L h-1]
    ##    Fair = u[7]        # Airflow inlet volumetric flowrate, [L h-1]
    ##    Tin = u[8]         # Temperature of bioreactor feed, [oC]
    ##    Tcin = u[9]        # Temperature of cooling agent inlet, [oC]

    # 1D Interpolation for Qin
    Qin = QinInterp(t)

    #Definition of model parameters
    #Kinetic parameters
    a1 = 0.05    # Ratkowsky parameter [oC-1 h-0.5]
    aP = 4.50    # Growth-associated parameter for ethanol production, [-]
    AP1 = 6.0    # Activation energy parameter for ethanol production, [oC]
    AP2 = 20.3   # Activation energy parameter for ethanol production, [oC]
    b1 = 0.035   # Parameter in the exponential expression of the maximum specific growth rate np.expression, [oC-1]
    b2 = 0.15    # Parameter in the exponential expression of the growth inhibitory ethanol concentration np.expression, [oC-1]
    b3 = 0.40    # Parameter in the exponential np.expression of the specific death rate expression,[oC-1]
    c1 = 0.38    # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
    c2 = 0.29    # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
    k1 = 3.00    # Parameter in the maximum specific growth rate expression, [oC]
    k2 = 55.0    # Parameter in the maximum specific growth rate expression, [oC]
    k3 = 60.0    # Parameter in the growth-inhibitory ethanol concentration expression, [oC]
    k4 = 50.0    # Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
    Pmaxb = 90   # Temperature-independent product inhibition constant, [g L-1]
    PmaxT = 90   # Maximum value of product inhibition constant due to temperature, [g L-1]
    Kdb = 0.025  # Basal specific cellular biomass death rate, [h-1]
    KdT = 30.00  # Maximum value of specific cellular biomass death rate due to temperature, [h-1]
    KSX = 5      # Glucose saturation constant for the specific growth rate, [g L-1]
    KOX = 0.0005 # Oxygen saturation constant for the specific growth rate, [g L-1]
    qOmax = 0.05 # Maximum specific oxygen consumption rate, [h-1]

    #Metabolic parameters
    YPS = 0.51 # Theoretical yield of ethanol on glucose, [gP gS-1]
    YXO = 0.97 # Theoretical yield of biomass on oxygen, [gX gO-1]
    YXS = 0.53 # Theoretical yield of biomass on glucose, [gX gS-1]

    #Physicochemical and thermodynamic parameters
    Chbr = 4.18     # Heat capacity of the mass of reaction, [J g-1 oC-1]
    Chc = 4.18      # Heat capacity of the cooling agent, [J g-1 oC-1]
    DeltaH = 518000 # Heat of reaction of fermentation, [J mol-1 O2]
    Tref = 20       # Reference temperature, [oC]
    KH = 200        # Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
    z = 0.792       # Oxygen compressibility factor, [-]
    R = 0.082       # Ideas gas constant, [L atm mol-1 oC-1]
    kla0 = 100      # Temperature-independent volumetric oxygen transfer coefficient, [h-1]
    KT = 360000     # Heat transfer coefficient, [J h-1 m-2 ??C-1]
    rho = 1080      # Density of the fermentation broth, [g L-1]
    rhoc = 1000     # Density of the cooling agent, [g L-1]
    MO = 32.0       # Molecular weight of oxygen (O2), [g mol-1]

    #Bioreactor design data
    AT = 1.0       # Bioreactor heat transfer area, [m2]
    V = 1800       # Bioreactor working volume, [L]
    Vcj = 50       # Cooling jacket volume, [L]
    Ogasin = 0.305 # Oxygen concentration in airflow inlet, [g L-1]

    #Definition of model variables
    #State variables
    Xt = x[0]      # Total cellular biomass, [g L-1]
    Xv = x[1]      # Viable cellular biomass, [g L-1]
    S = x[2]       # Substrate/Glucose concentration, [g L-1]
    P = x[3]       # Product/Ethanol concentration, [g L-1]
    Oliq = x[4]    # Dissolved oxygen concentration, [g L-1]
    Ogas = x[5]    # Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
    T = x[6]       # Temperature in the bioreactor, [oC]
    Tc = x[7]      # Temperature of the cooling agent in the jacket, [oC]
    Vl = x[8]      # Culture volume in the bioreactor, [L]
    Sf_cum = x[9]  # Cumulative amount of substrate/glucose fed to the bioreactor, [g]
    Time = x[10]   # Batch time, [h]

    # Definition of model equations
    # Kinetic rates
    # -----------------------------
    # Specific growth rate, [h-1] 
    mmax = ((a1*(T-k1))*(1-np.exp(b1*(T-k2))))**2
    Pmax = Pmaxb + PmaxT/(1-np.exp(-b2*(T-k3)))
    m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+np.exp(-(100-S)/1)) # Specific growth rate, [h-1] 
    if m1 >= 0:
        m = m1
    else:
        m=0.0
    # Non-growth-associated ethanol specific production rate, [h-1]
    if S > 0:
        bP = c1 * np.exp(-AP1/T) - c2 * np.exp(-AP2/T)  # Non-growth-associated ethanol specific production rate, [h-1]
    else:
        bP = 0.0
    qP = aP*m + bP
    # Ethanol consumption specific rate
    qS = m/YXS + qP/YPS
    # Oxygen consumption specific rate
    qO = qOmax*Oliq/YXO/(KOX + Oliq)
    # Specific biological deactivation rate of cell mass
    Kd = Kdb + KdT/(1+np.exp(-b3*(T-k4)))
    # Saturation concentration of oxygen in culture media
    Osat = z*Ogas*R*T/KH
    # Oxygen mass transfer coefficient
    kla = kla0*1.2**(T-20)
    # Volume of the gas phase in the bioreactor
    Vg = V - Vl

    #Material balances
    #-----------------
    # Volume of liquid culture
    dVl = Qin - Qe
    # Total cell mass
    dXt = m*Xv + Qin/Vl*(Xtin-Xt)
    # Total mass of biologically active cells
    dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv)
    # Glucose concentration
    dS = Qin/Vl*(Sin-S) - qS*Xv
    # Ethanol concentration
    dP = Qin/Vl*(-P) + qP*Xv
    # Disolved oxygen
    dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv
    # Oxygen gas phase
    dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg

    # Energy balances
    #---------------
    # Bioreactor temprature
    dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr
    # Cooling agent temperature
    dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc

    # Yields & Productivity
    #---------------------
    # Cumulative amount of glucose fed to the bioreactor
    dSf_cum = Sin*Qin
    dTime = 1

    # Definition of state derivatives vector
    # State derivatives
    dxdt = [dXt,dXv,dS,dP,dOliq,dOgas,dT,dTc,dVl,dSf_cum,dTime]

    return dxdt

T_val = np.linspace(7,40,21)
P_val = np.empty_like(T_val)
ux1 = ux0
for i,T in enumerate(T_val):
    ux1[9] = T
    x = odeint(ethanol,x0,tspan,args=tuple(ux1))
    P_val[i] = x[-3,3]

plt.figure(figsize=(6,3))
plt.plot(T_val, P_val, label='Maximize Product')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Tcin')
plt.tight_layout()
plt.savefig('bioreactor_Tcin.png',dpi=300)

F_val = np.linspace(10,1000,41)
P_val = np.empty_like(F_val)
ux2 = ux0
for i,F in enumerate(F_val):
    ux2[7] = F
    x = odeint(ethanol,x0,tspan,args=tuple(ux2))
    P_val[i] = x[-3,3]

plt.figure(figsize=(6,3))
plt.plot(F_val, P_val, label='Maximize Product')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Airflow [L/hr]')
plt.tight_layout()
plt.savefig('bioreactor_Fair.png',dpi=300)
plt.show()