У меня есть следующий скрипт, а переменная Q - это массив, вычисленный из другого скрипта, поэтому я хочу, чтобы этот Q менялся каждый раз, когда у меня решатель ODEINT: я хотел бы знать, нужно ли мне создавать отдельную функцию, которая содержит мой вектор Q или, если я должен ввести изменяющуюся время внутри функции rrn.
from scipy.integrate import odeint
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import opti_fedbatch_test_method2
# We import the parameters values from the following excel file
df = pd.read_excel("Kinetics_overview.xlsx")
# Here we define all the parameters used for the model (parameters have been extracted from literature)
# Glucose,Xylose,Biomass and Ethanol related parameters
nuMaxGlu = float(df['Value2'][df.index[df['Parameters'] == "nuMaxGlu"]]) # s-1
nuMaxXyl = float(df['Value2'][df.index[df['Parameters'] == "nuMaxXyl"]]) # s-1
Ks_Glu = float(df['Value2'][df.index[df['Parameters'] == "Ks_Glu"]]) # kg Glu m-3
Ks_Xyl = float(df['Value2'][df.index[df['Parameters'] == "Ks_Xyl"]]) # kg Xyl m-3
Ki_Glu = float(df['Value2'][df.index[df['Parameters'] == "Ki_Glu"]]) # kg Glu m-3
Ki_Xyl = float(df['Value2'][df.index[df['Parameters'] == "Ki_Xyl"]]) # kg Xyl m-3
Ki_GluXyl = float(df['Value2'][df.index[df['Parameters'] == "Ki_GluXyl"]]) # kg Glu m-3
Y_XGlu = float(df['Value2'][df.index[df['Parameters'] == "Y_XGlu"]]) # kg X/kg Glu
Y_XXyl = float(df['Value2'][df.index[df['Parameters'] == "Y_XXyl"]]) # kg X/kg Xyl
Ki_EtOHmaxGlu = float(df['Value2'][df.index[df['Parameters'] == "Ki_EtOHmaxGlu"]]) # kg Glu m-3
Ki_EtOHmaxXyl = float(df['Value2'][df.index[df['Parameters'] == "Ki_EtOHmaxXyl"]]) # kg Xyl m-3
Y_EtOHGlu = float(df['Value2'][df.index[df['Parameters'] == "Y_EtOHGlu"]]) # kg EtOH/kg Glu
Y_EtOHXyl = float(df['Value2'][df.index[df['Parameters'] == "Y_EtOHXyl"]]) # kg EtOH/kg Xyl
gammaG = float(df['Value2'][df.index[df['Parameters'] == "gammaG"]]) # no unit
gammaX = float(df['Value2'][df.index[df['Parameters'] == "gammaX"]]) # no unit
# Acetate-related parameters
nuHAcMax = float(df['Value2'][df.index[df['Parameters'] == "nuHAcMax"]]) # s-1
Ks_HAc = float(df['Value2'][df.index[df['Parameters'] == "Ks_HAc"]]) # kg HAc m-3
Ki_HAcGlu = float(df['Value2'][df.index[df['Parameters'] == "Ki_HAcGlu"]]) # kg HAc m-3
Ki_HAcXyl = float(df['Value2'][df.index[df['Parameters'] == "Ki_HAcXyl"]]) # kg HAc m-3
Y_HAcHMF = float(df['Value2'][df.index[df['Parameters'] == "Y_HAcHMF"]]) # kg Ac/kg HMF
# Furfural-related parameters
nuFurMax = float(df['Value2'][df.index[df['Parameters'] == "nuFurMax"]]) # s-1
Ks_Fur = float(df['Value2'][df.index[df['Parameters'] == "Ks_Fur"]]) # kg Furfural m-3
Ki_FurGlu = float(df['Value2'][df.index[df['Parameters'] == "Ki_FurGlu"]]) # kg Furfural m-3
Ki_FurXyl = float(df['Value2'][df.index[df['Parameters'] == "Ki_FurXyl"]]) # kg Furfural m-3
Ki_FurHMF = float(df['Value2'][df.index[df['Parameters'] == "Ki_FurHMF"]]) # kg Furfural m-3
Y_FurFA = float(df['Value2'][df.index[df['Parameters'] == "Y_FurFA"]]) # kg FA/kg Fur
# Furfuryl alcohol-related parameters
Ki_FAGlu = float(df['Value2'][df.index[df['Parameters'] == "Ki_FAGlu"]]) # kg FA m-3
Ki_FAXyl = float(df['Value2'][df.index[df['Parameters'] == "Ki_FAXyl"]]) # kg FA m-3
# HMF-related parameters
nuHMFMax = float(df['Value2'][df.index[df['Parameters'] == "nuHMFMax"]]) # s-1
Ks_HMF = float(df['Value2'][df.index[df['Parameters'] == "Ks_HMF"]]) # kg HMF m-3
Ki_HMFGlu = float(df['Value2'][df.index[df['Parameters'] == "Ki_HMFGlu"]]) # kg HMF m-3
Ki_HMFXyl = float(df['Value2'][df.index[df['Parameters'] == "Ki_HMFXyl"]]) # kg HMF m-3
# Initialization of the stoichiometric matrix
s = np.zeros((5, 8))
s[0, 0] = (Y_XGlu)
s[1, 0] = (Y_XXyl)
s[2, 0] = (0)
s[3, 0] = (0)
s[4, 0] = (0)
s[0, 1] = (-1)
s[1, 1] = (0)
s[2, 1] = (0)
s[3, 1] = (0)
s[4, 1] = (0)
s[0, 2] = (0)
s[1, 2] = (-1)
s[2, 2] = (0)
s[3, 2] = (0)
s[4, 2] = (0)
s[0, 3] = (Y_EtOHGlu)
s[1, 3] = (Y_EtOHXyl)
s[2, 3] = (0)
s[3, 3] = (0)
s[4, 3] = (0)
s[0, 4] = (0)
s[1, 4] = (0)
s[2, 4] = (-1)
s[3, 4] = (0)
s[4, 4] = (0)
s[0, 5] = (0)
s[1, 5] = (0)
s[2, 5] = (0)
s[3, 5] = (-1)
s[4, 5] = (Y_HAcHMF)
s[0, 6] = (0)
s[1, 6] = (0)
s[2, 6] = (0)
s[3, 6] = (0)
s[4, 6] = (-1)
s[0, 7] = (0)
s[1, 7] = (0)
s[2, 7] = (Y_FurFA)
s[3, 7] = (0)
s[4, 7] = (0)
# List of state variables
SV = ["Biomass", "Glucose", "Xylose", "Ethanol", "Furfural", "Acetic acid", "HMF", "Furfuryl alcohol", "Volume"]
# Initial values for the state variables
X0 = 10 # g/L
Glu0 = 20 # g/L
Xyl0 = 20 # g/L
EtOH0 = 0 # g/L
Fur0 = 1 # g/L
HAc0 = 1 # g/L
HMF0 = 1 # g/L
FA0 = 0 # g/L
V0 = 200 # L
# Optimal Q
Q = opti_fedbatch_test_method2.Q_opt_Xylrates[t] # L/h
# Initialization of the process vector used in the rxn function
rho = np.zeros(5)
# Initialization of the overall rates vector
r = np.zeros(8)
# Function that defines the equations of the model
def rxn(C, t):
# State variables inside of the Vector C which contains the concentration of the variables at
# every single time
X, Glu, Xyl, EtOH, Fur, HAc, HMF, FA, V = C
# Glucose uptake process
rho[0] = nuMaxGlu * X * (Glu / (Ks_Glu + Glu + ((Glu ** 2) / Ki_Glu)) *
(1 - (EtOH / Ki_EtOHmaxGlu) ** gammaG) *
(1 / (1 + (Fur / Ki_FurGlu))) *
(1 / (1 + (HAc / Ki_HAcGlu))) *
(1 / (1 + (HMF / Ki_HMFGlu))) *
(1 / (1 + (FA / Ki_FAGlu))))
# Xylose uptake process
rho[1] = nuMaxXyl * X * (Xyl / (Ks_Xyl + Xyl + ((Xyl ** 2) / Ki_Xyl)) *
(1 - (EtOH / Ki_EtOHmaxXyl) ** gammaX) *
(1 / (1 + (Fur / Ki_FurXyl))) *
(1 / (1 + (HAc / Ki_HAcXyl))) *
(1 / (1 + (HMF / Ki_HMFXyl))) *
(1 / (1 + (FA / Ki_FAXyl))) *
(1 / (1 + (Glu / Ki_GluXyl))))
# Fur uptake process
rho[2] = nuFurMax * X * (Fur / (Ks_Fur + Fur))
# HAc uptake process
rho[3] = nuHAcMax * X * (HAc / (Ks_HAc + HAc))
# HMF uptake process
rho[4] = nuHMFMax * X * (HMF / (Ks_HMF + HMF)) * (1 / (1 + (Fur / Ki_FurGlu)))
# Overall conversion rate (stoichiometric matrix * process vector)
r[0] = s[0, 0] * rho[0] + s[1, 0] * rho[1] + s[2, 0] * rho[2] + s[3, 0] * \
rho[3] + s[4, 0] * rho[4]
r[1] = s[0, 1] * rho[0] + s[1, 1] * rho[1] + s[2, 1] * rho[2] + s[3, 1] * \
rho[3] + s[4, 1] * rho[4]
r[2] = s[0, 2] * rho[0] + s[1, 2] * rho[1] + s[2, 2] * rho[2] + s[3, 2] * \
rho[3] + s[4, 2] * rho[4]
r[3] = s[0, 3] * rho[0] + s[1, 3] * rho[1] + s[2, 3] * rho[2] + s[3, 3] * \
rho[3] + s[4, 3] * rho[4]
r[4] = s[0, 4] * rho[0] + s[1, 4] * rho[1] + s[2, 4] * rho[2] + s[3, 4] * \
rho[3] + s[4, 4] * rho[4]
r[5] = s[0, 5] * rho[0] + s[1, 5] * rho[1] + s[2, 5] * rho[2] + s[3, 5] * \
rho[3] + s[4, 5] * rho[4]
r[6] = s[0, 6] * rho[0] + s[1, 6] * rho[1] + s[2, 6] * rho[2] + s[3, 6] * \
rho[3] + s[4, 6] * rho[4]
r[7] = s[0, 7] * rho[0] + s[1, 7] * rho[1] + s[2, 7] * rho[2] + s[3, 7] * \
rho[3] + s[4, 7] * rho[4]
# Mass balances for solving the ODE of each variable
dVdt = Q
dXdt = r[0] - X * (Q['a'](t) / V)
dGludt = r[1] + (Q / V) * (Glu0 - Glu)
dXyldt = r[2] + (Q / V) * (Xyl0 - Xyl)
dEtOHdt = r[3] - EtOH * (Q / V)
dFurdt = r[4] + (Q / V) * (Fur0 - Fur)
dHAcdt = r[5] + (Q / V) * (HAc0 - HAc)
dHMFdt = r[6] + (Q / V) * (Fur0 - Fur)
dFAdt = r[7] - FA * (Q / V)
return [dXdt, dGludt, dXyldt, dEtOHdt, dFurdt, dHAcdt, dHMFdt, dFAdt,
dVdt] # Here, we are solving the differential equation for each point.
# Here is performed the solving method
# Time vector to solve the ODE
t = np.linspace(0, 100, 101)
# Vector with the initial conditions of the state variables
C0 = [X0, Glu0, Xyl0, EtOH0, Fur0, HAc0, HMF0, FA0, V0]
# C vector containing the result of the solve
C = odeint(rxn, C0, t)
Concentrations = pd.DataFrame(C, columns=SV, index=t)
print("Initial Volume: " + str(V0) + " L")
print("Final Volume: " + str(C[-1, 8]) + " L")
print("Initial Xyl concentration: " + str(Xyl0) + " g/L")
print("Final Xyl concentration: " + str(C[-1, 2]) + " g/L")
print("Total amount of Xyl consumed: " + str((Xyl0 - C[-1, 2]) * C[-1, 8]) + " g")
print("Final EtOH amount: " + str(C[-1, 3] * C[-1, 8]) + " g")
print("Final EtOH concentration: " + str(C[-1, 3]) + " g/L")
# plotting code
plt.plot(t, C[:, 0], 'b-', label="X")
plt.plot(t, C[:, 1], 'g-', label="Glu")
plt.plot(t, C[:, 2], 'r-', label="Xyl")
plt.plot(t, C[:, 3], 'c-', label="EtOH")
plt.plot(t, C[:, 4], 'm-', label="Fur")
plt.plot(t, C[:, 5], 'y-', label="HAc")
plt.plot(t, C[:, 6], 'k-', label="HMF")
plt.plot(t, C[:, 7], '0.75', label="FA")
plt.xlabel("Time (s)")
plt.ylabel("Conc. (g/L)")
plt.legend(loc="upper right")
plt.show()
rho1 = np.zeros(100)
for i in range(100):
rho1[i] = nuMaxXyl * C[i, 0] * (C[i, 2] / (Ks_Xyl + C[i, 2] + ((C[i, 2] ** 2) / Ki_Xyl)) *
(1 - (C[i, 3] / Ki_EtOHmaxXyl) ** gammaX) *
(1 / (1 + (C[i, 4] / Ki_FurXyl))) *
(1 / (1 + (C[i, 5] / Ki_HAcXyl))) *
(1 / (1 + (C[i, 6] / Ki_HMFXyl))) *
(1 / (1 + (C[i, 7] / Ki_FAXyl))) *
(1 / (1 + (C[i, 1] / Ki_GluXyl))))