from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote=False)
# Create time vector: t=[0, 0.1, 0.2,...,36.9,37], [hours]
tm = np.linspace(0,37,371)
# Insert smaller time steps at the beginning
tm = np.insert(tm,1,[0.001,0.005,0.01,0.05])
m.time = tm
nt = len(tm)

# Define constants and parameters
#################################
# Kinetic Parameters
a1     = 0.05   # Ratkowsky parameter [oC-1 h-0.5]
aP     = 4.50   # Growth-associated parameter for EtOh production [-]
AP1    = 6.0    # Activation energy parameter for EtOh production [oC]
AP2    = 20.3   # Activation energy parameter for EtOh production [oC]
b1     = 0.035  # Parameter in the exponential expression of the maximum specific growth rate expression [oC-1]
b2     = 0.15   # Parameter in the exponential expression of the maximum specific growth rate expression [oC-1]
b3     = 0.40   # Parameter in the exponential expression of the specific death rate expression [oC-1]
c1     = 0.38   # Constant decoupling factor for EtOh [gP gX-1 h-1]
c2     = 0.29   # Constant decoupling factor for EtOh [gP gX-1 h-1]
k1     = 3      # Parameter in the maximum specific growth rate expression [oC]
k2     = 55     # Parameter in the maximum specific growth rate expression [oC]
k3     = 60     # Parameter in the growth-inhibitory EtOH concentration expression [oC]
k4     = 50     # 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     # 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 EtOH 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 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  # Ideal gas constant [L atm mol-1 oC-1]
kla0   = 100    # Temperature-independent volumetric oxygen transfer coefficient [-h]
KT     = 360000 # Heat transfer coefficient [J h-1 m-2 oC-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 [g mol-1]

# Bioreactor design data
AT     = 1      # 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]

# Define variables
##################
mi = m.Var(name='mi',lb=0)
# Qin is a step function:
#   Qin = Qin0 + 15H(t-5) + 5H(t-10) - 6H(t-20) - 14H(t-35)
#   where H(t-t0) heaviside function
Qin_step = np.zeros(nt)
Qin_step[np.where(tm>=5)]  += 15
Qin_step[np.where(tm>=10)] += 5
Qin_step[np.where(tm>=20)] -= 6
Qin_step[np.where(tm>=35)] -= 14
Qin  = m.Param(value=Qin_step, name='Qin')
# Fixed variables, they are constant throughout the time horizon
Xtin = m.FV(value=0, name='Xtin')
Xvin = m.FV(value=0, name='Xvin')
Qe   = m.FV(value=0, name='Qe')
Sin  = m.FV(value=400, lb=0, ub=1500)
Pin  = m.FV(value=0, name='Pin')
Fc   = m.FV(value=40, name='Fc')
Fair = m.FV(value=60000, name='Fair')
Tin  = m.FV(value=30, name='Tin')
Tcin = m.FV(value=15, name='Tcin')

Vl     = m.Var(value=1000, name='Vl') # lb=-0.0, ub=0.75*V
Xt     = m.Var(value=0.1,  name='Xt') # lb=-0.0, ub=10
Xv     = m.Var(value=0.1,  name='Xv') # lb=-0.0, ub=10
S      = m.Var(value=50,  name='S')   # lb=+0.0, ub=10000
P      = m.Var(value=0, name='P')
Ol     = m.Var(value=0.0065, name= 'Ol')
Og     = m.Var(value=0.305, name='Og')
T      = m.Var(value=30, name='T')    # lb=20, ub=40
Tc     = m.Var(value=20, name='Tc')   # lb=0, ub=30
Sf_cum = m.Var(value=0, name='Sf_cum')

# Define algebraic equations
############################
# Specific growth rate of cell mass
mimax = m.Intermediate(((a1*(T-k1))*(1-m.exp(b1*(T-k2))))** 2)
Pmax = m.Intermediate(Pmaxb + PmaxT/(1-m.exp(-b2*(T-k3))))
m.Equation(mi == mimax * (S / (KSX+S)) * (Ol/(KOX + Ol)) \
                 * (1 - P/Pmax) * (1 / (1+m.exp(-(100-S)))))
# Specific production rate of EtOH
bP = m.Intermediate(c1*m.exp(-AP1/T) - c2*m.exp(-AP2/T))
qP = m.Intermediate(aP*mi + bP)
# Specific consumption rate of glucose
qS = m.Intermediate(mi/YXS + qP/YPS)
# Specific consumption rate of oxygen
qO = m.Intermediate(qOmax*Ol/YXO/(KOX+Ol))
# Specific biological deactivation rate of cell mass
Kd = m.Intermediate(Kdb + KdT/(1+m.exp(-b3*(T-k4))))
# Saturation concentration of oxygen in culture media
Ostar = m.Intermediate(z*Og*R*T/KH)
# Oxygen mass transfer coefficient
kla = m.Intermediate(kla0*1.2**(T-20))
# Bioreactor phases equation
Vg = m.Intermediate(V - Vl)

# Define differential equations
###############################
m.Equation(Vl.dt()     == Qin - Qe)
m.Equation(Vl*Xt.dt()  == Qin*(Xtin-Xt) + mi*Vl*Xv)
m.Equation(Vl*Xv.dt()  == Qin*(Xvin-Xv) + Xv*Vl*(mi-Kd))
m.Equation(Vl*S.dt()   == Qin*(Sin-S) - qS*Vl*Xv)
m.Equation(Vl*P.dt()   == Qin*(Pin - P) + qP*Vl*Xv)
m.Equation(Vl*Ol.dt()  == Qin*(Ostar-Ol)   + Vl*kla*(Ostar-Ol) - qO*Vl*Xv)
m.Equation(Vg*Og.dt()  == Fair*(Ogasin-Og) - Vl*kla*(Ostar-Ol) + Og*(Qin-Qe))
m.Equation(Vl*T.dt()   == Qin*(Tin-T) - Tref*(Qin-Qe) \
                          + Vl*qO*Xv*deltaH/MO/rho/Chbr - KT*AT*(T-Tc)/rho/Chbr)
m.Equation(Vcj*Tc.dt() == Fc*(Tcin - Tc) + KT*AT*(T-Tc)/rhoc/Chc)
m.Equation(Sf_cum.dt() == Qin*Sin)

# solve ODE
m.options.SOLVER= 1
m.options.IMODE = 7
m.options.NODES = 3
# m.open_folder()
m.solve(disp=False)

# Plot results
plt.figure(1)
plt.title('Total & Viable Cellular Biomass')
plt.plot(m.time, Xv.value, label='Xv')
plt.plot(m.time, Xt.value, label='Xt')
plt.legend()
plt.ylabel('Biomass concentration [g/L]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()

plt.figure(2)
plt.title('Substrate (S) & Product (P) concentration')
plt.plot(m.time, S.value, label='S')
plt.plot(m.time, P.value, label='P')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Time [h]')
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()

plt.figure(3)
plt.title('Bioreactor & Cooling jacket temperature')
plt.plot(m.time, T.value, label='T')
plt.plot(m.time, Tc.value, label='Tc')
plt.legend()
plt.ylabel('Temperature [oC]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()

fig4, ax = plt.subplots()
ax.title.set_text('Dissolved & Gaseous Oxygen concentration')
lns1 = ax.plot(m.time, Ol.value, label='[Oliq]', color='c')
ax.set_xlabel('Time [h]')
ax.set_ylabel('Oliq [g/L]', color='c')
ax.minorticks_on()
ax2 = ax.twinx()
lns2 = ax2.plot(m.time, Og.value, label='[Ogas]', color='y')
ax2.set_ylabel('Ogas [g/L]', color='y')
ax2.minorticks_on()
lns = lns1 + lns2
labs = [l.get_label() for l in lns]
ax.legend(lns, labs, loc='best')
ax.grid()
fig4.tight_layout()
plt.figure(4)

plt.figure(5)
plt.title('Feeding Policy')
plt.plot(m.time, Qin.value, label='Qin')
plt.legend()
plt.ylabel('Qin [L/h]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()

plt.figure(6)
plt.title('Check >=0 Constraints')
plt.plot(tm,bP.value,label='bP')
plt.plot(tm,mi.value,label='mi')
plt.legend(); plt.grid()

plt.show()