我有一个化学反应系统，其中速率方程代表间歇反应器模型。该模型是使用 SciPysolve_ivp函数求解的 ODE 系统。我的示例（见下文）目前包含 5 个反应，但我需要从整体方案中添加更多反应。使用我目前的方法，跟踪化学物质、质量分数、速率参数等是很乏味的。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Parameters
# -----------------------------------------------------------------------------

species = (
    'CELL', 'CELLA', 'CH2OHCHO', 'CHOCHO', 'CH3CHO', 'C6H6O3', 'C2H5CHO',
    'CH3OH', 'CH2O', 'CO', 'GCO', 'CO2', 'H2', 'H2O', 'GCH2O', 'HCOOH',
    'CH2OHCH2CHO', 'CH4', 'GH2', 'CHAR', 'C6H10O5', 'GCOH2', 'GMSW', 'HCE1',
    'HCE2'
)

# time vector to evaluate reaction rates [s]
time = np.linspace(0, 2.0)

# Solve batch reactor system of equations
# -----------------------------------------------------------------------------

def dcdt_debiagi(t, y):
    """
    Rate equations for Debiagi 2018 biomass pyrolysis kinetics.
    """
    R = 1.9859          # universal gas constant [cal/(mol K)]
    T = 773.15          # temperature [K]

    cell = y[0]         # cellulose mass fraction [-]
    cella = y[1]        # active cellulose mass fraction [-]
    gmsw = y[22]        # softwood hemicellulose mass fraction [-]

    # Cellulose reactions and rate constants
    # 1) CELL -> CELLA
    # 2) CELLA -> 0.40 CH2OHCHO + 0.03 CHOCHO + 0.17 CH3CHO + 0.25 C6H6O3 + 0.35 C2H5CHO + 0.20 CH3OH + 0.15 CH2O + 0.49 CO + 0.05 GCO + 0.43 CO2 + 0.13 H2 + 0.93 H2O + 0.05 GCH2O + 0.02 HCOOH + 0.05 CH2OHCH2CHO + 0.05 CH4 + 0.1 GH2 + 0.66 CHAR
    # 3) CELLA -> C6H10O5
    # 4) CELL -> 4.45 H2O + 5.45 CHAR + 0.12 GCOH2 + 0.18 GCH2O + 0.25 GCO + 0.125 GH2 + 0.125 H2
    k1 = 1.5e14 * np.exp(-47_000 / (R * T))
    k2 = 2.5e6 * np.exp(-19_100 / (R * T))
    k3 = 3.3 * T * np.exp(-10_000 / (R * T))
    k4 = 9e7 * np.exp(-31_000 / (R * T))

    # Hemicellulose reactions and rate constants
    # 5) GMSW -> 0.70 HCE1 + 0.30 HCE2
    k5 = 1e10 * np.exp(-31_000 / (R * T))

    # species reaction rate equations where r = dc/dt
    # mass fractions associated with each species are also given
    r_CELL = -k1 * cell
    r_CELLA = k1 * cell - k2 * cella - k3 * cella
    r_CH2OHCHO = k2 * cella * 0.1481
    r_CHOCHO = k2 * cella * 0.0107
    r_CH3CHO = k2 * cella * 0.0462
    r_C6H6O3 = k2 * cella * 0.1944
    r_C2H5CHO = k2 * cella * 0.1254
    r_CH3OH = k2 * cella * 0.0395
    r_CH2O = k2 * cella * 0.0278
    r_CO = k2 * cella * 0.0846
    r_GCO = k2 * cella * 0.008637 + k4 * cell * 0.04319
    r_CO2 = k2 * cella * 0.1167
    r_H2 = k2 * cella * 0.001616 + k4 * cell * 0.001554
    r_H2O = k2 * cella * 0.1033 + k4 * cell * 0.4944
    r_GCH2O = k2 * cella * 0.009259 + k4 * cell * 0.03333
    r_HCOOH = k2 * cella * 0.005677
    r_CH2OHCH2CHO = k2 * cella * 0.02284
    r_CH4 = k2 * cella * 0.004947
    r_GH2 = k2 * cella * 0.001243 + k4 * cell * 0.001554
    r_CHAR = k2 * cella * 0.04889 + k4 * cell * 0.4037
    r_C6H10O5 = k3 * cella
    r_GCOH2 = k4 * cell * 0.02222
    r_GMSW = -k5 * gmsw
    r_HCE1 = k5 * gmsw * 0.7
    r_HCE2 = k5 * gmsw * 0.3

    # system of ODEs
    dcdt = (
        r_CELL, r_CELLA, r_CH2OHCHO, r_CHOCHO, r_CH3CHO, r_C6H6O3, r_C2H5CHO,
        r_CH3OH, r_CH2O, r_CO, r_GCO, r_CO2, r_H2, r_H2O, r_GCH2O, r_HCOOH,
        r_CH2OHCH2CHO, r_CH4, r_GH2, r_CHAR, r_C6H10O5, r_GCOH2, r_GMSW, r_HCE1,
        r_HCE2
    )
    return dcdt

# initial mass fractions [-]
y0 = np.zeros(len(species))
y0[0] = 0.6
y0[22] = 0.4

# solution for system of equations for a batch reactor
sol = solve_ivp(dcdt_debiagi, (time[0], time[-1]), y0, t_eval=time)

# Print
# ----------------------------------------------------------------------------

print('--- Initial ---')
print(f'CELL        {sol.y[0][0]}')
print(f'GMSW        {sol.y[22][0]}')

print('--- Final ---')
for i, sp in enumerate(species):
    print(f'{sp:11} {sol.y[i][-1]:.4f}')

理想情况下，我想定义一个数组中的所有反应，解析该数组，并输出一个矩阵，SciPy 可以使用该矩阵solve_ivp来确定每个时间步的物种浓度。我正在考虑如下所示的内容：

R = 1.9859      # universal gas constant [cal/(mol K)]
T = 773.15      # temperature [K]

reactions = [
    'CELL -> CELLA',
    'CELLA -> 0.40 CH2OHCHO + 0.03 CHOCHO + 0.17 CH3CHO + 0.25 C6H6O3 + 0.35 C2H5CHO + 0.20 CH3OH + 0.15 CH2O + 0.49 CO + 0.05 GCO + 0.43 CO2 + 0.13 H2 + 0.93 H2O + 0.05 GCH2O + 0.02 HCOOH + 0.05 CH2OHCH2CHO + 0.05 CH4 + 0.1 GH2 + 0.66 CHAR',
    'CELLA -> C6H10O5'
]

names = {
    'CELL': 'C6H10O5',
    'CELLA': 'C6H10O5',
    'GCO': 'CO',
    'GCH2O': 'CH2O',
    'GH2': 'H2',
    'CHAR': 'C'
}

arrhenius = [
    (1.5e14, 0, 47_000),
    (2.5e6, 0, 19_100),
    (3.3, 1, 10_000)
]

如何定义一个化学方程式系统，solve_ivp无需为每种化学物质手动编写每个速率方程式即可实现？