I'm trying to perform fitting for bioprocess model consisting of 5 ODEs system and 6 algebraic equations using symfit library. The main issue is that I have experimental data for variables F and kLa which are used in ODEs equations but are not the state variables, and can't figure out how to implement this variables into symfit.
EDIT:system of equations
Below is my code:
import numpy as np
from symfit import parameters, variables, Fit, D, ODEModel
import matplotlib.pyplot as plt
import pandas as pd
columns=['t','OD','S','cond','DO','V','F','kLa']
Dane = pd.read_csv('exp_data.txt',sep='\t', names=columns)
#Experimental data
timepoints=np.array(Dane['t']) # time points
X_data=np.array(Dane['OD']) # X measurements
S_data=np.array(Dane['S']) # S measurements
A_data=np.array(Dane['cond']) # A measurements
DO_data=np.array(Dane['DO']) # O measurements
V_data=np.array(Dane['V']) # V measurements
F_data=np.array(Dane['F']) # F measurements
kLa_data=np.array(Dane['kLa']) # kLa measurements
#initial values
t0=Dane['t'][0]
X0=Dane['OD'][0]
S0=Dane['S'][0]
A0=Dane['cond'][0]
DO0=Dane['DO'][0]
V0=Dane['V'][0]
F0=Dane['F'][0]
kLa0=Dane['kLa'][0]
Sf=500
# Define variables
t, X, S, A, DO, V, F, kLa = variables('t X S A DO V F kLa')
# Define parameters with bounds
Kap, Ksa, Ko, Ks, Kia, Kis, p_Amax, q_Amax, qm, q_Smax, Yoa, Yxa, Yem, Yos, Yxsof = parameters(
'Kap Ksa Ko Ks Kia Kis p_Amax q_Amax qm q_Smax Yoa Yxa Yem Yos Yxsof',
bounds={
'Kap': (0, None), # Lower bound of 0 for Kap, no upper bound
'Ksa': (0, None),
'Ko': (0, None),
'Ks': (0, None),
'Kia': (0, None),
'Kis': (0, None),
'p_Amax': (0, None),
'q_Amax': (0, None),
'qm': (0, None),
'q_Smax': (0, None),
'Yoa': (0, None),
'Yxa': (0, None),
'Yem': (0, None),
'Yos': (0, None),
'Yxsof': (0, None),
}
)
# Define state variables and parameters for the algebraic equations
qs = (q_Smax * S / (S + Ks)) * Kia / (Kia + A)
qsof = (p_Amax * qs / (qs + Kap))
qsox = (qs - qsof) * DO / (DO + Ko)
qsa = (q_Amax * A / (A + Ksa)) * (Kis / (qs + Kis))
qo = (qsox - qm) * Yos + qsa * Yoa
u = (qsox - qm) * Yem + qsof * Yxsof + qsa * Yxa
# Define the system of differential equations
dx1 = u * X - F * X / V
dx2 = (F * (Sf - S) / V) - qs * X
dx3 = qsa * X - (F * A / V)
dx4 = kLa * (100 - DO) - qo * X * 700
dx5 = F
# Define the ODE model
model={
D(X, t): dx1,
D(S, t): dx2,
D(A, t): dx3,
D(DO, t): dx4,
D(V, t): dx5,
}
ode_model = ODEModel(model,initial={t:t0, X:X0, S:S0, A:A0, DO:DO0, V:V0})
# Fit the model to the data
fit = Fit(ode_model, t=timepoints, X=X_data, S=S_data, A=A_data, DO=DO_data, V=V_data)
result = fit.execute()
# Display the result
print(result)
# Access the fitted parameters
fitted_params = result.params
print(fitted_params)
Running this code gives the following Traceback and Error:
Traceback (most recent call last):
File "D:\Python praca\Optymalizacaj\dopasowanie parametrow hodowli.py", line 91, in <module>
result = fit.execute()
File "C:\ProgramData\Anaconda3\lib\site-packages\symfit\core\fit.py", line 573, in execute
minimizer_ans = self.minimizer.execute(**minimize_options)
File "C:\ProgramData\Anaconda3\lib\site-packages\symfit\core\minimizers.py", line 400, in execute
return super(ScipyGradientMinimize, self).execute(jacobian=jacobian, **minimize_options)
File "C:\ProgramData\Anaconda3\lib\site-packages\symfit\core\minimizers.py", line 339, in execute
ans = minimize(
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\_minimize.py", line 612, in minimize
return _minimize_bfgs(fun, x0, args, jac, callback, **options)
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\optimize.py", line 1101, in _minimize_bfgs
sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\optimize.py", line 261, in _prepare_scalar_function
sf = ScalarFunction(fun, x0, args, grad, hess,
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 76, in __init__
self._update_fun()
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 166, in _update_fun
self._update_fun_impl()
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 73, in update_fun
self.f = fun_wrapped(self.x)
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 70, in fun_wrapped
return fun(x, *args)
File "C:\ProgramData\Anaconda3\lib\site-packages\symfit\core\objectives.py", line 308, in __call__
evaluated_func = super(LeastSquares, self).__call__(
File "C:\ProgramData\Anaconda3\lib\site-packages\symfit\core\objectives.py", line 93, in __call__
result = self.model(**key2str(parameters))._asdict()
File "C:\ProgramData\Anaconda3\lib\site-packages\symfit\core\models.py", line 1256, in __call__
return ModelOutput(self.keys(), self.eval_components(*args, **kwargs))
File "C:\ProgramData\Anaconda3\lib\site-packages\symfit\core\models.py", line 1194, in eval_components
ans_bigger = solve_ivp(
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 576, in solve_ivp
message = solver.step()
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 181, in step
success, message = self._step_impl()
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\lsoda.py", line 148, in _step_impl
solver._y, solver.t = integrator.run(
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ode.py", line 1346, in run
y1, t, istate = self.runner(*args)
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 138, in fun
return self.fun_single(t, y)
File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 20, in fun_wrapped
return np.asarray(fun(t, y), dtype=dtype)
File "C:\ProgramData\Anaconda3\lib\site-packages\numpy\core\_asarray.py", line 85, in asarray
return array(a, dtype, copy=False, order=order)
File "C:\ProgramData\Anaconda3\lib\site-packages\sympy\core\expr.py", line 327, in __float__
raise TypeError("can't convert expression to float")
TypeError: can't convert expression to float
As I understand, probably variables F and kLa can't be converted to the numeric values, but I don't know how to overcome this. I have tried adding initial conditions for these variables in ode_model as also below code with no success:
fit = Fit(ode_model, t=timepoints, X=X_data, S=S_data, A=A_data, DO=DO_data, V=V_data, F=F_data, kLa=kLa_data)
How to pass this variables to fitting in a proper way? Appreciate any suggestions and comments.
EDIT: test data, the same format as in my source file:
0 0.9 2 18.4 99 6250 18.9 121.956
1 1.3 1.2 18.59 64.7 6268.9 12.6 121.843
2 2.5 1.2 18.92 44.7 6281.5 26.6 122.069
3 4.4 1.0 19.53 45.2 6308.1 47.7 138.115
4 7.8 1.1 19.93 37.7 6355.8 88.2 150.356
5 12.8 1.1 20.73 36 6444.0 116.8 164.471
6 21 1.0 22.15 37.6 6560.8 188.6 176.819
6.5 26 0.9 23.26 42.5 6655.1 139.2 185.219
7 29 0.9 24.54 37.8 6724.7 73.2 187.883
8 35 0 22.83 57.8 6797.9 60.8 192.626
9 44 0 21.01 46.4 6858.7 74.0 192.626
10 52 0 18.76 44 6932.7 75.6 192.678
11 58 0 17.2 49.8 7008.3 67.2 192.739
12 65 0 15.67 66.9 7075.5 53.2 192.709
13 72 0 14.63 83 7128.7 39.2 192.709
14 75 0 14.09 89.8 7167.9 36.4 192.791
14.5 76 0 13.93 92.3 7186.1 25.2 192.804
15 76 0 13.84 92.9 7198.7 25.2 192.722