如何使用 symfit 将非状态变量放入 ODE 系统优化模型中？

Question

我正在尝试使用 symfit 库对由 5 个 ODE 系统和 6 个代数方程组成的生物过程模型进行拟合。主要问题是我有变量 F 和 kLa 的实验数据，它们在 ODE 方程中使用，但不是状态变量，并且无法弄清楚如何将这些变量实现到 symfit 中。

编辑：方程组

下面是我的代码：

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)

运行此代码会出现以下回溯和错误：

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

据我了解，变量 F 和 kLa 可能无法转换为数值，但我不知道如何克服这个问题。我尝试在 ode_model 中为这些变量添加初始条件，如下代码所示，但没有成功：

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)

如何以正确的方式将此变量传递给拟合？感谢任何建议和意见。

编辑：测试数据，与我的源文件中的格式相同：

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

您需要拟合 F 和 kLa，然后根据以这种方式获得的参数创建符号表达式。

我相应地修改了代码，现在它可以工作了：

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)

输出：

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)])

拟合 F 的图：

拟合 kLa 的图：

请注意，由于某种原因，命名变量 S 会导致 symfit 崩溃。我重命名了变量 Z，一切正常。