How to put variable which is not a state variable in a ODEs system optimisation model using symfit?

Question

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

Answer 1

You need to fit F and kLa and then create symbolic expressions out of the parameters obtained that way.

I modified the code accordingly and it now works:

import numpy as np
from symfit import parameters, variables, Fit, D, ODEModel
import symfit as sf
import matplotlib.pyplot as plt
import pandas as pd
from functools import reduce
from scipy.optimize import curve_fit

columns=['t','OD','S','cond','DO','V','F','kLa']
Dane = pd.read_csv('exp_data.csv', 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

# Fit F as a sum of sine functions
def f(t, *args, sin=np.sin):
    terms = map(lambda a, b, c: a * sin(b * t + c), *[args[i*3: (i+1)*3] for i in range(len(args) // 3)])
    result = reduce(lambda a, b: a + b, terms) + args[-1]
    return result

p0 = np.array([-20.84503322, -33.7246875 ,  43.35807928,   1.64346027, \
                0.99209514,   0.53846615,  51.60185929,  -0.31651934, \
              -20.76987162,  75.26997873])

fit_result_F = curve_fit(f, timepoints, F_data, p0)
t = np.linspace(0, 15, 500)
plt.figure()
plt.scatter(timepoints, F_data)
p_F = fit_result_F[0]
plt.plot(t, f(t, *p_F))

# Fit kLa to a sigmoid like function
def sigmoid(t, a, b, c, exp=np.exp):
    return a/(1+exp(-t+b)) + c

p0 = (70, 1, 120)
fit_result_kLa = curve_fit(sigmoid, timepoints, kLa_data, p0)
t = np.linspace(0, 15, 500)
plt.figure()
plt.scatter(timepoints, kLa_data)
p_kLa = fit_result_kLa[0]
plt.plot(t, sigmoid(t, *p_kLa))

#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 Z 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

F = f(t, *p_F, sin=sf.sin)
kLa = sigmoid(t, *p_kLa, exp=sf.exp)

# Define the ODE model
# Define the system of differential equations
model={
    D(X, t): u * X - F * X / V,
    D(S, t): (F * (Sf - S) / V) - qs * X,
    D(A, t): qsa * X - (F * A / V),
    D(DO, t): kLa * (100 - DO) - qo * X * 700,
    D(V, t): F,
}

ode_model = ODEModel(model, {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, Z=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)

Output:

Parameter Value        Standard Deviation
Kap       1.000001e+00 7.939781e-01
Kia       1.000002e+00 5.158846e-01
Kis       1.000000e+00 3.633269e-01
Ko        9.999968e-01 1.013783e+00
Ks        9.999972e-01 1.283476e+00
Ksa       1.000008e+00 4.712817e-01
Yem       9.999952e-01 9.719080e-01
Yoa       9.999919e-01 4.922445e-01
Yos       9.999952e-01 5.265178e-01
Yxa       1.000006e+00 4.020917e-01
Yxsof     9.999960e-01 1.096456e+00
p_Amax    9.999929e-01 4.132026e-01
q_Amax    9.999962e-01 3.848603e-01
q_Smax    9.999978e-01 8.411599e-01
qm        9.999990e-01 8.769199e-01
Status message         Desired error not necessarily achieved due to precision loss.
Number of iterations   1
Objective              <symfit.core.objectives.LeastSquares object at 0x7f7a1affcf10>
Minimizer              <symfit.core.minimizers.BFGS object at 0x7f7a1b024d90>

Goodness of fit qualifiers:
chi_squared            149758.6523249151
objective_value        74879.32616245755
r_squared              0.9299590820921964
OrderedDict([('Kap', 1.000000929381065), ('Kia', 1.0000015696787052), ('Kis', 1.0000004109657716), ('Ko', 0.9999967819815748), ('Ks', 0.9999971637966988), ('Ksa', 1.000008181072036), ('Yem', 0.999995216882535), ('Yoa', 0.999991885855043), ('Yos', 0.9999952459083982), ('Yxa', 1.0000055360671372), ('Yxsof', 0.9999960317243858), ('p_Amax', 0.999992898633446), ('q_Amax', 0.9999962302286587), ('q_Smax', 0.9999978034646714), ('qm', 0.9999989848526278)])

Plot of fitted F:

Plot of fitted kLa:

Note that for some reason, naming the variable S causes symfit to crash. I renamed the variable Z and all is OK.