Core ADM1 Model

pyadm1.core.adm1.ADM1

ADM1da – 41-state ODE system.

Implements the sub-fraction disintegration/hydrolysis approach, temperature-dependent kinetics, and modified inhibition kinetics described in Schlattmann (2011); SIMBA# biogas 4.2 Tutorial.

Usage

from pyadm1 import Feedstock fs = Feedstock([...], feeding_freq=24) adm = ADM1(fs, V_liq=1200, V_gas=216, T_ad=315.15) adm.set_influent_dataframe(fs.get_influent_dataframe(Q=[11.4, 6.1])) adm.create_influent(Q=[11.4, 6.1], i=0) state0 = [0.01] * STATE_SIZE dydt = adm.ADM_ODE(0.0, state0)

Source code in pyadm1/core/adm1.py

class ADM1:
    """
    ADM1da – 41-state ODE system.

    Implements the sub-fraction disintegration/hydrolysis approach,
    temperature-dependent kinetics, and modified inhibition kinetics
    described in Schlattmann (2011); SIMBA# biogas 4.2 Tutorial.

    Usage
    -----
    >>> from pyadm1 import Feedstock
    >>> fs = Feedstock([...], feeding_freq=24)
    >>> adm = ADM1(fs, V_liq=1200, V_gas=216, T_ad=315.15)
    >>> adm.set_influent_dataframe(fs.get_influent_dataframe(Q=[11.4, 6.1]))
    >>> adm.create_influent(Q=[11.4, 6.1], i=0)
    >>> state0 = [0.01] * STATE_SIZE
    >>> dydt = adm.ADM_ODE(0.0, state0)
    """

    def __init__(
        self,
        feedstock,
        V_liq: float = 1977.0,
        V_gas: float = 304.0,
        T_ad: float = 308.15,
    ) -> None:
        """
        Initialize the ADM1 model.

        Parameters
        ----------
        feedstock : Feedstock
            Feedstock object (used by ``create_influent`` when no external
            influent DataFrame has been provided via ``set_influent_dataframe``).
        V_liq : float
            Liquid digester volume [m³].
        V_gas : float
            Gas headspace volume [m³].
        T_ad : float
            Operating temperature [K] (default 308.15 K = 35 °C).
        """
        # --- Reactor volumes ---
        self.V_liq = V_liq
        self._V_gas = V_gas
        self._V_ad = V_liq + V_gas

        # --- Temperature ---
        self._T_ad = T_ad

        # --- Physical constants ---
        self._R = 0.08314  # bar·m³·kmol⁻¹·K⁻¹
        self._T_base = 298.15  # 25 °C reference
        self._p_atm = 1.013

        self._RT = self._R * self._T_ad
        self._p_ext = self._p_atm - 0.0084147 * np.exp(0.054 * (self._T_ad - 273.15))

        # --- Feedstock / influent ---
        self._feedstock = feedstock
        self._Q: Optional[List[float]] = None
        self._state_input: Optional[List[float]] = None

        # --- Calibration overrides ---
        self._calibration_params: dict = {}

        # --- Result-tracking lists ---
        self._Q_GAS: List[float] = []
        self._Q_CH4: List[float] = []
        self._Q_CO2: List[float] = []
        self._Q_H2O: List[float] = []
        self._P_GAS: List[float] = []
        self._pH_l: List[float] = []
        self._FOSTAC: List[float] = []
        self._AcvsPro: List[float] = []
        self._VFA: List[float] = []
        self._TAC: List[float] = []

        # --- Temperature-corrected kinetics (reused across ODE calls) ---
        base_kinetic = ADMParams.get_kinetic_params()
        theta = ADMParams.get_temperature_factors()
        self._kinetic = ADMParams.apply_temperature_corrections(base_kinetic, theta, T_ad)
        # Snapshot of the temperature-corrected defaults so calibration
        # overrides can be reverted cleanly by clear_calibration_parameters().
        self._kinetic_default = dict(self._kinetic)

        # --- Stoichiometric / fraction / inhibition parameters ---
        self._stoich = ADMParams.get_stoichiometric_params()
        self._fractions = ADMParams.get_product_fractions()
        self._inhib_params = ADMParams.get_inhibition_params(self._R, self._T_base, T_ad)

        # --- Gas parameters ---
        p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2 = ADMParams.getADMgasparams(self._R, self._T_base, T_ad)
        self._p_gas_h2o = p_gas_h2o
        self._k_p = k_p
        self._k_L_a = k_L_a
        self._K_H_co2 = K_H_co2
        self._K_H_ch4 = K_H_ch4
        self._K_H_h2 = K_H_h2
        self._K_H_co2_default = K_H_co2
        self._K_H_ch4_default = K_H_ch4
        self._K_H_h2_default = K_H_h2

        # Gas volume reference conditions (theta = 20 °C, 1 atm) — ADM1da convention
        self._T_gas_norm = 293.15
        self._P_gas_norm = 1.01325
        self._NQ = 1000.0 * self._P_gas_norm / (self._R * self._T_gas_norm)

        # Pressure threshold for linearised gas-outlet pipe formula [bar].
        self._dP_min_pipe = 1.0e-3

        # Optional external influent DataFrame.
        self._influent_df: Optional[pd.DataFrame] = None

        # Influent / sludge densities (used by ``create_influent`` only).
        self._rho_in: float = 1000.0
        self._rho_sludge: float = 1000.0

        # When ``_Q_out_override`` is not None, ADM_ODE uses it for the sludge
        # outflow instead of the default ``q_in - q_S_loss`` formula. The
        # caller is then responsible for advancing V_liq between ODE steps.
        # ``_q_S_loss_last`` caches the most recent q_S_loss so the caller
        # can integrate the volume balance externally.
        self._Q_out_override: Optional[float] = None
        self._q_S_loss_last: float = 0.0

    # ------------------------------------------------------------------
    # Public read-only properties
    # ------------------------------------------------------------------

    @property
    def model_name(self) -> str:
        """Short identifier for this model variant."""
        return "ADM1"

    @property
    def T_ad(self) -> float:
        """Operating temperature [K]."""
        return self._T_ad

    @property
    def feedstock(self):
        """Feedstock object."""
        return self._feedstock

    @property
    def Q_GAS(self) -> List[float]:
        """History of total biogas flow rate [m^3/d]."""
        return self._Q_GAS

    @property
    def Q_CH4(self) -> List[float]:
        """History of methane flow rate [m^3/d]."""
        return self._Q_CH4

    @property
    def Q_CO2(self) -> List[float]:
        """History of carbon dioxide flow rate [m^3/d]."""
        return self._Q_CO2

    @property
    def Q_H2O(self) -> List[float]:
        """History of water-vapour flow rate [m^3/d]."""
        return self._Q_H2O

    @property
    def P_GAS(self) -> List[float]:
        """History of total headspace pressure [bar]."""
        return self._P_GAS

    @property
    def pH_l(self) -> List[float]:
        """History of liquid-phase pH."""
        return self._pH_l

    @property
    def VFA_TA(self) -> List[float]:
        """History of the VFA/TA (FOS/TAC) ratio."""
        return self._FOSTAC

    @property
    def AcvsPro(self) -> List[float]:
        """History of the acetate-to-propionate ratio."""
        return self._AcvsPro

    @property
    def VFA(self) -> List[float]:
        """History of total volatile fatty acid concentration [kg HAc-eq/m^3]."""
        return self._VFA

    @property
    def TAC(self) -> List[float]:
        """History of total alkalinity (TAC)."""
        return self._TAC

    def get_state_size(self) -> int:
        """Return the number of state variables: 41."""
        return STATE_SIZE

    # ------------------------------------------------------------------
    # Calibration parameter API
    # ------------------------------------------------------------------

    def set_calibration_parameters(self, parameters: dict) -> None:
        """Set calibration overrides for kinetic / gas parameters.

        Kinetic-rate keys (``k_m_*``, ``k_hyd_*``, ``k_dis_*``, ``Y_*``,
        ``K_S_*``, ``k_dec_*``) are applied directly into the live kinetic
        dict consumed by the ODE. Henry-constant keys (``K_H_co2``,
        ``K_H_ch4``, ``K_H_h2``) overwrite the corresponding attributes.
        ``k_p`` and ``k_L_a`` continue to be looked up from
        ``self._calibration_params`` at their respective call sites.
        """
        self._calibration_params.update(parameters)
        for key, value in parameters.items():
            if key in self._kinetic_default:
                self._kinetic[key] = float(value)
            elif key == "K_H_co2":
                self._K_H_co2 = float(value)
            elif key == "K_H_ch4":
                self._K_H_ch4 = float(value)
            elif key == "K_H_h2":
                self._K_H_h2 = float(value)

    def clear_calibration_parameters(self) -> None:
        """Clear all calibration parameters and restore defaults."""
        self._calibration_params = {}
        self._kinetic = dict(self._kinetic_default)
        self._K_H_co2 = self._K_H_co2_default
        self._K_H_ch4 = self._K_H_ch4_default
        self._K_H_h2 = self._K_H_h2_default

    def get_calibration_parameters(self) -> dict:
        """Return a copy of the current calibration parameters."""
        return self._calibration_params.copy()

    # ------------------------------------------------------------------
    # Gas calculation
    # ------------------------------------------------------------------

    def calc_gas(
        self,
        pi_Sh2: float,
        pi_Sch4: float,
        pi_Sco2: float,
        pTOTAL: float,
    ) -> Tuple[float, float, float, float, float]:
        """
        Calculate biogas production rates, normalised to standard reference
        conditions (theta = 20 °C, P = 1.01325 bar; ADM1da convention).

        Returns
        -------
        q_gas, q_ch4, q_co2, q_h2o, p_gas
        """
        k_p = self._k_p
        if self._calibration_params.get("k_p") is not None:
            k_p = float(self._calibration_params["k_p"])

        p_total_wet = pTOTAL + self._p_gas_h2o
        q_gas = max(
            k_p * (p_total_wet - self._p_ext) / (self._RT / 1000.0 * self._NQ) * self.V_liq,
            0.0,
        )

        p_gas = pi_Sh2 + pi_Sch4 + pi_Sco2
        p_gas_wet = p_gas + self._p_gas_h2o
        if p_gas_wet > 0.0:
            q_ch4 = max(q_gas * (pi_Sch4 / p_gas_wet), 0.0)
            q_co2 = max(q_gas * (pi_Sco2 / p_gas_wet), 0.0)
            q_h2o = max(q_gas * (self._p_gas_h2o / p_gas_wet), 0.0)
        else:
            q_ch4 = 0.0
            q_co2 = 0.0
            q_h2o = 0.0

        return q_gas, q_ch4, q_co2, q_h2o, p_gas

    # ------------------------------------------------------------------
    # Influent setup
    # ------------------------------------------------------------------

    def set_influent_dataframe(self, df: pd.DataFrame) -> None:
        """
        Store an external ADM1 influent DataFrame.

        The DataFrame must have columns matching ``INFLUENT_COLUMNS`` (37
        liquid-state columns + Q).  Once set, ``create_influent()`` reads
        from this DataFrame instead of deriving values from the feedstock.
        """
        missing = [c for c in INFLUENT_COLUMNS if c not in df.columns]
        if missing:
            raise ValueError(f"Influent DataFrame missing columns: {missing}")
        self._influent_df = df.reset_index(drop=True)

    def set_influent_density(self, rho_in: float, rho_sludge: float = 1000.0) -> None:
        """Retained for API compatibility; stored but no longer affects Q_out."""
        self._rho_in = float(rho_in)
        self._rho_sludge = float(rho_sludge)

    def create_influent(
        self,
        Q: List[float],
        i: int,
        rho: Optional[List[float]] = None,
    ) -> None:
        """
        Build the influent vector for time step *i*.

        If an external influent DataFrame has been set via
        ``set_influent_dataframe()``, that DataFrame is used.  Otherwise
        the feedstock object is used to derive the influent.
        """
        if hasattr(self._feedstock, "actual_Q"):
            Q_actual = self._feedstock.actual_Q(Q)
        else:
            Q_actual = list(Q)
        self._Q = Q_actual
        if rho is not None and len(rho) == len(Q_actual):
            q_total = sum(Q_actual)
            if q_total > 0.0:
                self._rho_in = sum(q * r for q, r in zip(Q_actual, rho)) / q_total

        if self._influent_df is not None:
            max_i = len(self._influent_df) - 1
            row = self._influent_df.iloc[min(i, max_i)]
            self._state_input = row[INFLUENT_COLUMNS[:-1]].tolist()
        else:
            df = self._feedstock.get_influent_dataframe(Q)
            max_i = len(df) - 1
            row = df.iloc[min(i, max_i)]
            self._state_input = row[INFLUENT_COLUMNS[:-1]].tolist()

    # ------------------------------------------------------------------
    # ODE system
    # ------------------------------------------------------------------

    def ADM_ODE(self, t: float, state: List[float]) -> Tuple[float, ...]:
        """
        Compute dy/dt for the 41-element ADM1 state vector.

        The system is autonomous; *t* is present only to satisfy the scipy
        solver interface.
        """
        # ---- Unpack state ------------------------------------------------
        S_su = state[_IDX_S_SU]
        S_aa = state[_IDX_S_AA]
        S_fa = state[_IDX_S_FA]
        S_va = state[_IDX_S_VA]
        S_bu = state[_IDX_S_BU]
        S_pro = state[_IDX_S_PRO]
        S_ac = state[_IDX_S_AC]
        S_h2 = state[_IDX_S_H2]
        S_ch4 = state[_IDX_S_CH4]
        S_co2 = state[_IDX_S_CO2]
        S_nh4 = state[_IDX_S_NH4]
        S_I = state[_IDX_S_I]

        X_PS_ch = state[_IDX_X_PS_CH]
        X_PS_pr = state[_IDX_X_PS_PR]
        X_PS_li = state[_IDX_X_PS_LI]
        X_PF_ch = state[_IDX_X_PF_CH]
        X_PF_pr = state[_IDX_X_PF_PR]
        X_PF_li = state[_IDX_X_PF_LI]
        X_S_ch = state[_IDX_X_S_CH]
        X_S_pr = state[_IDX_X_S_PR]
        X_S_li = state[_IDX_X_S_LI]
        X_I = state[_IDX_X_I]

        X_su = state[_IDX_X_SU]
        X_aa = state[_IDX_X_AA]
        X_fa = state[_IDX_X_FA]
        X_c4 = state[_IDX_X_C4]
        X_pro = state[_IDX_X_PRO]
        X_ac = state[_IDX_X_AC]
        X_h2 = state[_IDX_X_H2]

        S_cation = state[_IDX_S_CATION]
        S_anion = state[_IDX_S_ANION]
        S_va_ion = state[_IDX_S_VA_ION]
        S_bu_ion = state[_IDX_S_BU_ION]
        S_pro_ion = state[_IDX_S_PRO_ION]
        S_ac_ion = state[_IDX_S_AC_ION]
        S_hco3 = state[_IDX_S_HCO3]
        S_nh3 = state[_IDX_S_NH3]

        p_gas_h2 = state[_IDX_P_H2]
        p_gas_ch4 = state[_IDX_P_CH4]
        p_gas_co2 = state[_IDX_P_CO2]
        pTOTAL = state[_IDX_P_TOTAL]

        # ---- Influent and hydraulic flow ---------------------------------
        q_ad = float(np.sum(self._Q)) if self._Q is not None else 0.0
        s_in = self._state_input if self._state_input is not None else [0.0] * 37

        # ---- Parameters --------------------------------------------------
        k = self._kinetic
        st = self._stoich
        fr = self._fractions
        ip = self._inhib_params

        f_fa_li = st["f_fa_li"]
        f_ch_bac = st["f_ch_bac"]
        f_pr_bac = st["f_pr_bac"]
        f_li_bac = st["f_li_bac"]
        f_p_bac = st["f_p_bac"]
        fXI_PS = st["fXI_PS"]
        fXI_PF = st["fXI_PF"]
        fSI = st["fSI_hyd"]

        Y_su = k["Y_su"]
        Y_aa = k["Y_aa"]
        Y_fa = k["Y_fa"]
        Y_c4 = k["Y_c4"]
        Y_pro = k["Y_pro"]
        Y_ac = k["Y_ac"]
        Y_h2 = k["Y_h2"]

        # ---- pH (Newton–Raphson charge balance) --------------------------
        S_H = self._calc_ph(
            S_nh4,
            S_nh3,
            S_hco3,
            S_ac_ion,
            S_pro_ion,
            S_bu_ion,
            S_va_ion,
            S_cation,
            S_anion,
            ip["K_w"],
        )
        S_H = max(S_H, 1.0e-14)

        # ---- Inhibition factors ------------------------------------------
        S_IN = S_nh4 + S_nh3
        I_IN = S_IN / (ip["K_S_IN"] + S_IN + 1.0e-20)

        I_pH_aa = self._pH_inhib(S_H, ip["K_pH_aa"], n=1)
        I_pH_fa = self._pH_inhib(S_H, ip["K_pH_aa"], n=2)
        I_pH_c4 = self._pH_inhib(S_H, ip["K_pH_aa"], n=2)
        I_pH_pro = self._pH_inhib(S_H, ip["K_pH_aa"], n=2)
        I_pH_ac = self._pH_inhib(S_H, ip["K_pH_ac"], n=3)
        I_pH_h2 = self._pH_inhib(S_H, ip["K_pH_h2"], n=3)

        I_h2_fa = ip["K_I_h2_fa"] / (ip["K_I_h2_fa"] + S_h2)
        I_h2_c4 = ip["K_I_h2_c4"] / (ip["K_I_h2_c4"] + S_h2)
        I_h2_pro = ip["K_I_h2_pro"] / (ip["K_I_h2_pro"] + S_h2)

        K_nh3_ac_sq = ip["K_I_nh3"] ** 2
        K_nh3_pro_sq = ip["K_I_nh3_pro"] ** 2
        S_nh3_sq = S_nh3 * S_nh3
        I_nh3 = K_nh3_ac_sq / (K_nh3_ac_sq + S_nh3_sq)
        I_nh3_pro = K_nh3_pro_sq / (K_nh3_pro_sq + S_nh3_sq)

        K_co2_h2_sq = ip["K_S_co2_h2"] * ip["K_S_co2_h2"]
        S_co2_sq = S_co2 * S_co2
        I_co2_h2 = S_co2_sq / (K_co2_h2_sq + S_co2_sq + 1.0e-30)

        Ka_pro = ip["K_a_pro"]
        S_HPr = (S_pro / 112.0) * S_H / (S_H + Ka_pro + 1.0e-20)
        I_HPr = ip["K_IH_pro"] / (ip["K_IH_pro"] + S_HPr + 1.0e-20)

        Ka_ac = ip["K_a_ac"]
        S_HAc = (S_ac / 64.0) * S_H / (S_H + Ka_ac + 1.0e-20)
        I_HAc = ip["K_IH_ac"] / (ip["K_IH_ac"] + S_HAc + 1.0e-20)

        I_su = I_pH_aa * I_IN
        I_aa = I_pH_aa * I_IN
        I_fa = I_pH_fa * I_IN * I_h2_fa
        I_c4 = I_pH_c4 * I_IN * I_h2_c4
        I_pro = I_pH_pro * I_IN * I_h2_pro * I_HPr * I_nh3_pro
        I_ac = I_pH_ac * I_IN * I_nh3 * I_HAc
        I_h2 = I_pH_h2 * I_IN * I_co2_h2

        # ---- Process rates -----------------------------------------------
        Rho_dis_PS_ch = k["k_dis_PS"] * X_PS_ch
        Rho_dis_PS_pr = k["k_dis_PS"] * X_PS_pr
        Rho_dis_PS_li = k["k_dis_PS"] * X_PS_li
        Rho_dis_PF_ch = k["k_dis_PF"] * X_PF_ch
        Rho_dis_PF_pr = k["k_dis_PF"] * X_PF_pr
        Rho_dis_PF_li = k["k_dis_PF"] * X_PF_li

        Rho_hyd_ch = k["k_hyd_ch"] * X_S_ch
        Rho_hyd_pr = k["k_hyd_pr"] * X_S_pr
        Rho_hyd_li = k["k_hyd_li"] * X_S_li

        Rho_su = k["k_m_su"] * S_su / (k["K_S_su"] + S_su + 1.0e-20) * X_su * I_su
        Rho_aa = k["k_m_aa"] * S_aa / (k["K_S_aa"] + S_aa + 1.0e-20) * X_aa * I_aa
        Rho_fa = k["k_m_fa"] * S_fa / (k["K_S_fa"] + S_fa + 1.0e-20) * X_fa * I_fa

        S_vbu = S_va + S_bu + 1.0e-20
        Rho_c4_va = k["k_m_c4"] * S_va / (k["K_S_c4"] + S_va + 1.0e-20) * X_c4 * (S_va / S_vbu) * I_c4
        Rho_c4_bu = k["k_m_c4"] * S_bu / (k["K_S_c4"] + S_bu + 1.0e-20) * X_c4 * (S_bu / S_vbu) * I_c4

        Rho_pro = k["k_m_pro"] * S_pro / (k["K_S_pro"] + S_pro + 1.0e-20) * X_pro * I_pro
        Rho_ac = k["k_m_ac"] * S_ac / (k["K_S_ac"] + S_ac + 1.0e-20) * X_ac * I_ac
        Rho_h2 = k["k_m_h2"] * S_h2 / (k["K_S_h2"] + S_h2 + 1.0e-20) * X_h2 * I_h2

        Rho_dec_su = k["k_dec_su"] * X_su
        Rho_dec_aa = k["k_dec_aa"] * X_aa
        Rho_dec_fa = k["k_dec_fa"] * X_fa
        Rho_dec_c4 = k["k_dec_c4"] * X_c4
        Rho_dec_pro = k["k_dec_pro"] * X_pro
        Rho_dec_ac = k["k_dec_ac"] * X_ac
        Rho_dec_h2 = k["k_dec_h2"] * X_h2
        sum_decay = Rho_dec_su + Rho_dec_aa + Rho_dec_fa + Rho_dec_c4 + Rho_dec_pro + Rho_dec_ac + Rho_dec_h2

        # ---- Acid-base rates --------------------------------------------
        k_AB = ip["k_A_B"]
        Rho_A_va = k_AB * (S_va_ion * S_H - ip["K_a_va"] * (S_va - S_va_ion))
        Rho_A_bu = k_AB * (S_bu_ion * S_H - ip["K_a_bu"] * (S_bu - S_bu_ion))
        Rho_A_pro = k_AB * (S_pro_ion * S_H - ip["K_a_pro"] * (S_pro - S_pro_ion))
        Rho_A_ac = k_AB * (S_ac_ion * S_H - ip["K_a_ac"] * (S_ac - S_ac_ion))
        Rho_A_co2 = k_AB * (S_hco3 * S_H - ip["K_a_co2"] * (S_co2 - S_hco3))
        Rho_A_IN = k_AB * (S_nh3 * S_H - ip["K_a_IN"] * (S_nh4 - S_nh3))

        # ---- Gas transfer rates -----------------------------------------
        k_L_a = self._k_L_a
        if self._calibration_params.get("k_L_a") is not None:
            k_L_a = float(self._calibration_params["k_L_a"])

        Rho_T_h2 = k_L_a * (S_h2 - 16.0 * p_gas_h2 / (self._RT * self._K_H_h2)) * (self.V_liq / self._V_gas)
        Rho_T_ch4 = k_L_a * (S_ch4 - 64.0 * p_gas_ch4 / (self._RT * self._K_H_ch4)) * (self.V_liq / self._V_gas)
        S_co2_free = max(S_co2 - S_hco3, 0.0)
        Rho_T_co2 = k_L_a * (S_co2_free - p_gas_co2 / (self._RT * self._K_H_co2)) * (self.V_liq / self._V_gas)

        k_p = self._k_p
        if self._calibration_params.get("k_p") is not None:
            k_p = float(self._calibration_params["k_p"])
        Rho_T_11 = max(
            k_p * (pTOTAL + self._p_gas_h2o - self._p_ext) * (self.V_liq / self._V_gas),
            0.0,
        )

        # ---- Carbon stoichiometry coefficients for S_co2 balance --------
        C = st
        s_hyd_ch = -C["C_ch"] + (1.0 - fSI) * C["C_su"] + fSI * C["C_I_s"]
        s_hyd_pr = -C["C_pr"] + (1.0 - fSI) * C["C_aa"] + fSI * C["C_I_s"]
        s_hyd_li = -C["C_li"] + (1.0 - fSI) * (f_fa_li * C["C_fa"] + (1.0 - f_fa_li) * C["C_su"]) + fSI * C["C_I_s"]

        s_su = (
            -C["C_su"]
            + (1.0 - Y_su) * (fr["f_bu_su"] * C["C_bu"] + fr["f_pro_su"] * C["C_pro"] + fr["f_ac_su"] * C["C_ac"])
            + Y_su * C["C_bac"]
        )
        s_aa = (
            -C["C_aa"]
            + (1.0 - Y_aa)
            * (fr["f_va_aa"] * C["C_va"] + fr["f_bu_aa"] * C["C_bu"] + fr["f_pro_aa"] * C["C_pro"] + fr["f_ac_aa"] * C["C_ac"])
            + Y_aa * C["C_bac"]
        )
        s_fa = -C["C_fa"] + (1.0 - Y_fa) * 0.7 * C["C_ac"] + Y_fa * C["C_bac"]
        s_c4_va = -C["C_va"] + (1.0 - Y_c4) * (0.54 * C["C_pro"] + 0.31 * C["C_ac"]) + Y_c4 * C["C_bac"]
        s_c4_bu = -C["C_bu"] + (1.0 - Y_c4) * 0.8 * C["C_ac"] + Y_c4 * C["C_bac"]
        s_pro = -C["C_pro"] + (1.0 - Y_pro) * 0.57 * C["C_ac"] + Y_pro * C["C_bac"]
        s_ac = -C["C_ac"] + (1.0 - Y_ac) * C["C_ch4"] + Y_ac * C["C_bac"]
        s_h2 = (1.0 - Y_h2) * C["C_ch4"] + Y_h2 * C["C_bac"]
        s_dec = -C["C_bac"] + f_ch_bac * C["C_ch"] + f_pr_bac * C["C_pr"] + f_li_bac * C["C_li"] + f_p_bac * C["C_I_x"]

        Sigma = (
            s_hyd_ch * Rho_hyd_ch
            + s_hyd_pr * Rho_hyd_pr
            + s_hyd_li * Rho_hyd_li
            + s_su * Rho_su
            + s_aa * Rho_aa
            + s_fa * Rho_fa
            + s_c4_va * Rho_c4_va
            + s_c4_bu * Rho_c4_bu
            + s_pro * Rho_pro
            + s_ac * Rho_ac
            + s_h2 * Rho_h2
            + s_dec * sum_decay
        )

        # ---- Differential equations -------------------------------------
        D_in = q_ad / self.V_liq

        # Sludge volume loss (hydrolysis-rate volume balance, after Schlattmann 2011).
        _q_S_loss = self.V_liq * (
            Rho_hyd_ch * (0.9375 / 1550.0) + Rho_hyd_pr * (0.6125 / 1370.0) + Rho_hyd_li * (0.3474 / 920.0)
        )
        self._q_S_loss_last = float(_q_S_loss)

        # If the caller has supplied an outflow (e.g. via a level controller
        # that owns the volume balance), use it. Otherwise enforce volume
        # conservation by setting Q_out = Q_in - q_S_loss.
        if self._Q_out_override is not None:
            _Q_out = max(float(self._Q_out_override), 0.0)
        else:
            _Q_out = max(q_ad - _q_S_loss, 0.0)
        D_out = _Q_out / self.V_liq

        # --- Dissolved (0–11) ---
        diff_S_su = (
            D_in * s_in[0] - D_out * S_su + (1.0 - fSI) * Rho_hyd_ch + (1.0 - fSI) * (1.0 - f_fa_li) * Rho_hyd_li - Rho_su
        )
        diff_S_aa = D_in * s_in[1] - D_out * S_aa + (1.0 - fSI) * Rho_hyd_pr - Rho_aa
        diff_S_fa = D_in * s_in[2] - D_out * S_fa + (1.0 - fSI) * f_fa_li * Rho_hyd_li - Rho_fa
        diff_S_va = D_in * s_in[3] - D_out * S_va + (1.0 - Y_aa) * fr["f_va_aa"] * Rho_aa - Rho_c4_va
        diff_S_bu = (
            D_in * s_in[4]
            - D_out * S_bu
            + (1.0 - Y_su) * fr["f_bu_su"] * Rho_su
            + (1.0 - Y_aa) * fr["f_bu_aa"] * Rho_aa
            - Rho_c4_bu
        )
        diff_S_pro = (
            D_in * s_in[5]
            - D_out * S_pro
            + (1.0 - Y_su) * fr["f_pro_su"] * Rho_su
            + (1.0 - Y_aa) * fr["f_pro_aa"] * Rho_aa
            + (1.0 - Y_c4) * 0.54 * Rho_c4_va
            - Rho_pro
        )
        diff_S_ac = (
            D_in * s_in[6]
            - D_out * S_ac
            + (1.0 - Y_su) * fr["f_ac_su"] * Rho_su
            + (1.0 - Y_aa) * fr["f_ac_aa"] * Rho_aa
            + (1.0 - Y_fa) * 0.7 * Rho_fa
            + (1.0 - Y_c4) * 0.31 * Rho_c4_va
            + (1.0 - Y_c4) * 0.8 * Rho_c4_bu
            + (1.0 - Y_pro) * 0.57 * Rho_pro
            - Rho_ac
        )
        diff_S_h2 = (
            D_in * s_in[7]
            - D_out * S_h2
            + (1.0 - Y_su) * fr["f_h2_su"] * Rho_su
            + (1.0 - Y_aa) * fr["f_h2_aa"] * Rho_aa
            + (1.0 - Y_fa) * 0.3 * Rho_fa
            + (1.0 - Y_c4) * 0.15 * Rho_c4_va
            + (1.0 - Y_c4) * 0.2 * Rho_c4_bu
            + (1.0 - Y_pro) * 0.43 * Rho_pro
            - Rho_h2
            - self._V_gas / self.V_liq * Rho_T_h2
        )
        diff_S_ch4 = (
            D_in * s_in[8]
            - D_out * S_ch4
            + (1.0 - Y_ac) * Rho_ac
            + (1.0 - Y_h2) * Rho_h2
            - self._V_gas / self.V_liq * Rho_T_ch4
        )
        diff_S_co2 = D_in * s_in[9] - D_out * S_co2 - Sigma - self._V_gas / self.V_liq * Rho_T_co2 + Rho_A_co2

        N_bac = st["N_bac"]
        N_aa = st["N_aa"]
        N_I = st["N_I"]
        diff_S_nh4 = (
            D_in * s_in[10]
            - D_out * S_nh4
            - Y_su * N_bac * Rho_su
            + (N_aa - Y_aa * N_bac) * Rho_aa
            - Y_fa * N_bac * Rho_fa
            - Y_c4 * N_bac * Rho_c4_va
            - Y_c4 * N_bac * Rho_c4_bu
            - Y_pro * N_bac * Rho_pro
            - Y_ac * N_bac * Rho_ac
            - Y_h2 * N_bac * Rho_h2
            + (N_bac - f_pr_bac * N_aa - f_p_bac * N_I) * sum_decay
            + Rho_A_IN
        )
        diff_S_I = D_in * s_in[11] - D_out * S_I + fSI * (Rho_hyd_ch + Rho_hyd_pr + Rho_hyd_li)

        # --- Particulate sub-fractions (12–21) ---
        diff_X_PS_ch = D_in * s_in[12] - D_out * X_PS_ch - Rho_dis_PS_ch
        diff_X_PS_pr = D_in * s_in[13] - D_out * X_PS_pr - Rho_dis_PS_pr
        diff_X_PS_li = D_in * s_in[14] - D_out * X_PS_li - Rho_dis_PS_li
        diff_X_PF_ch = D_in * s_in[15] - D_out * X_PF_ch - Rho_dis_PF_ch
        diff_X_PF_pr = D_in * s_in[16] - D_out * X_PF_pr - Rho_dis_PF_pr
        diff_X_PF_li = D_in * s_in[17] - D_out * X_PF_li - Rho_dis_PF_li
        diff_X_S_ch = (
            D_in * s_in[18]
            - D_out * X_S_ch
            + (1.0 - fXI_PS) * Rho_dis_PS_ch
            + (1.0 - fXI_PF) * Rho_dis_PF_ch
            + f_ch_bac * sum_decay
            - Rho_hyd_ch
        )
        diff_X_S_pr = (
            D_in * s_in[19]
            - D_out * X_S_pr
            + (1.0 - fXI_PS) * Rho_dis_PS_pr
            + (1.0 - fXI_PF) * Rho_dis_PF_pr
            + f_pr_bac * sum_decay
            - Rho_hyd_pr
        )
        diff_X_S_li = (
            D_in * s_in[20]
            - D_out * X_S_li
            + (1.0 - fXI_PS) * Rho_dis_PS_li
            + (1.0 - fXI_PF) * Rho_dis_PF_li
            + f_li_bac * sum_decay
            - Rho_hyd_li
        )
        diff_X_I = (
            D_in * s_in[21]
            - D_out * X_I
            + fXI_PS * (Rho_dis_PS_ch + Rho_dis_PS_pr + Rho_dis_PS_li)
            + fXI_PF * (Rho_dis_PF_ch + Rho_dis_PF_pr + Rho_dis_PF_li)
            + f_p_bac * sum_decay
        )

        # --- Biomass (22–28) ---
        diff_X_su = D_in * s_in[22] - D_out * X_su + Y_su * Rho_su - Rho_dec_su
        diff_X_aa = D_in * s_in[23] - D_out * X_aa + Y_aa * Rho_aa - Rho_dec_aa
        diff_X_fa = D_in * s_in[24] - D_out * X_fa + Y_fa * Rho_fa - Rho_dec_fa
        diff_X_c4 = D_in * s_in[25] - D_out * X_c4 + Y_c4 * Rho_c4_va + Y_c4 * Rho_c4_bu - Rho_dec_c4
        diff_X_pro = D_in * s_in[26] - D_out * X_pro + Y_pro * Rho_pro - Rho_dec_pro
        diff_X_ac = D_in * s_in[27] - D_out * X_ac + Y_ac * Rho_ac - Rho_dec_ac
        diff_X_h2 = D_in * s_in[28] - D_out * X_h2 + Y_h2 * Rho_h2 - Rho_dec_h2

        # --- Charge balance (29–36) ---
        diff_S_cation = D_in * s_in[29] - D_out * S_cation
        diff_S_anion = D_in * s_in[30] - D_out * S_anion
        diff_S_va_ion = D_in * s_in[31] - D_out * S_va_ion - Rho_A_va
        diff_S_bu_ion = D_in * s_in[32] - D_out * S_bu_ion - Rho_A_bu
        diff_S_pro_ion = D_in * s_in[33] - D_out * S_pro_ion - Rho_A_pro
        diff_S_ac_ion = D_in * s_in[34] - D_out * S_ac_ion - Rho_A_ac
        diff_S_hco3 = D_in * s_in[35] - D_out * S_hco3 - Rho_A_co2
        diff_S_nh3 = D_in * s_in[36] - D_out * S_nh3 - Rho_A_IN

        # --- Gas phase (37–40) ---
        diff_p_h2 = Rho_T_h2 * self._RT / 16.0 - p_gas_h2 / pTOTAL * Rho_T_11
        diff_p_ch4 = Rho_T_ch4 * self._RT / 64.0 - p_gas_ch4 / pTOTAL * Rho_T_11
        diff_p_co2 = Rho_T_co2 * self._RT - p_gas_co2 / pTOTAL * Rho_T_11
        diff_pTOT = self._RT / 16.0 * Rho_T_h2 + self._RT / 64.0 * Rho_T_ch4 + self._RT * Rho_T_co2 - Rho_T_11

        return (
            diff_S_su,
            diff_S_aa,
            diff_S_fa,
            diff_S_va,
            diff_S_bu,
            diff_S_pro,
            diff_S_ac,
            diff_S_h2,
            diff_S_ch4,
            diff_S_co2,
            diff_S_nh4,
            diff_S_I,
            diff_X_PS_ch,
            diff_X_PS_pr,
            diff_X_PS_li,
            diff_X_PF_ch,
            diff_X_PF_pr,
            diff_X_PF_li,
            diff_X_S_ch,
            diff_X_S_pr,
            diff_X_S_li,
            diff_X_I,
            diff_X_su,
            diff_X_aa,
            diff_X_fa,
            diff_X_c4,
            diff_X_pro,
            diff_X_ac,
            diff_X_h2,
            diff_S_cation,
            diff_S_anion,
            diff_S_va_ion,
            diff_S_bu_ion,
            diff_S_pro_ion,
            diff_S_ac_ion,
            diff_S_hco3,
            diff_S_nh3,
            diff_p_h2,
            diff_p_ch4,
            diff_p_co2,
            diff_pTOT,
        )

    # ------------------------------------------------------------------
    # Result-tracking helper
    # ------------------------------------------------------------------

    def print_params_at_current_state(self, state: List[float]) -> None:
        """Compute and store process indicators (pH, gas) from the current state."""
        ip = self._inhib_params
        S_H = self._calc_ph(
            state[_IDX_S_NH4],
            state[_IDX_S_NH3],
            state[_IDX_S_HCO3],
            state[_IDX_S_AC_ION],
            state[_IDX_S_PRO_ION],
            state[_IDX_S_BU_ION],
            state[_IDX_S_VA_ION],
            state[_IDX_S_CATION],
            state[_IDX_S_ANION],
            ip["K_w"],
        )
        pH = -np.log10(max(S_H, 1.0e-14))
        self._track_pH(round(pH, 1))

        q_gas, q_ch4, q_co2, q_h2o, p_gas = self.calc_gas(
            state[_IDX_P_H2],
            state[_IDX_P_CH4],
            state[_IDX_P_CO2],
            state[_IDX_P_TOTAL],
        )
        self._track_gas(q_gas, q_ch4, q_co2, q_h2o, p_gas)

    def resume_from_broken_simulation(self, Q_CH4: List[float]) -> None:
        """Re-populate the methane tracking list after a simulation restart."""
        for q in Q_CH4:
            self._Q_CH4.append(q)

    def _track_pH(self, pH: float) -> None:
        """Append pH to the history; pad on the very first call to keep length >= 2."""
        self._pH_l.append(pH)
        if len(self._pH_l) < 2:
            self._pH_l.append(self._pH_l[-1])

    def _track_gas(self, q_gas, q_ch4, q_co2, q_h2o, p_gas) -> None:
        """Append gas flows and headspace pressure to their history lists; pad on the first call."""
        self._Q_GAS.append(q_gas)
        self._Q_CH4.append(q_ch4)
        self._Q_CO2.append(q_co2)
        self._Q_H2O.append(q_h2o)
        self._P_GAS.append(p_gas)
        if len(self._Q_GAS) < 2:
            for _ in range(2):
                self._Q_GAS.append(q_gas)
                self._Q_CH4.append(q_ch4)
                self._Q_CO2.append(q_co2)
                self._Q_H2O.append(q_h2o)
                self._P_GAS.append(p_gas)

    # ------------------------------------------------------------------
    # Private helpers
    # ------------------------------------------------------------------

    @staticmethod
    def _pH_inhib(S_H: float, K_pH: float, n: int = 1) -> float:
        """Hill-type pH inhibition factor: I = K_pH^n / (K_pH^n + S_H^n)."""
        Kn = K_pH**n
        SHn = S_H**n
        return Kn / (Kn + SHn)

    @staticmethod
    def _calc_ph(
        S_nh4: float,
        S_nh3: float,
        S_hco3: float,
        S_ac_ion: float,
        S_pro_ion: float,
        S_bu_ion: float,
        S_va_ion: float,
        S_cation: float,
        S_anion: float,
        K_w: float,
        max_iter: int = 50,
    ) -> float:
        """Solve the charge balance for [H+] using Newton–Raphson iteration."""
        vfa_anions = S_ac_ion / 64.0 + S_pro_ion / 112.0 + S_bu_ion / 160.0 + S_va_ion / 208.0
        fixed = S_cation - S_anion + (S_nh4 - S_nh3) - S_hco3 - vfa_anions

        S_H = 1.0e-7
        for _ in range(max_iter):
            f = fixed + S_H - K_w / (S_H + 1.0e-30)
            df = 1.0 + K_w / (S_H + 1.0e-30) ** 2
            delta = -f / df
            S_H = max(1.0e-14, S_H + delta)
            if abs(delta) < 1.0e-15:
                break
        return S_H

Attributes

AcvsPro property

History of the acetate-to-propionate ratio.

P_GAS property

History of total headspace pressure [bar].

Q_CH4 property

History of methane flow rate [m^3/d].

Q_CO2 property

History of carbon dioxide flow rate [m^3/d].

Q_GAS property

History of total biogas flow rate [m^3/d].

Q_H2O property

History of water-vapour flow rate [m^3/d].

TAC property

History of total alkalinity (TAC).

T_ad property

Operating temperature [K].

VFA property

History of total volatile fatty acid concentration [kg HAc-eq/m^3].

VFA_TA property

History of the VFA/TA (FOS/TAC) ratio.

feedstock property

Feedstock object.

model_name property

Short identifier for this model variant.

pH_l property

History of liquid-phase pH.

Functions

ADM_ODE(t, state)

Compute dy/dt for the 41-element ADM1 state vector.

The system is autonomous; t is present only to satisfy the scipy solver interface.

Source code in pyadm1/core/adm1.py

def ADM_ODE(self, t: float, state: List[float]) -> Tuple[float, ...]:
    """
    Compute dy/dt for the 41-element ADM1 state vector.

    The system is autonomous; *t* is present only to satisfy the scipy
    solver interface.
    """
    # ---- Unpack state ------------------------------------------------
    S_su = state[_IDX_S_SU]
    S_aa = state[_IDX_S_AA]
    S_fa = state[_IDX_S_FA]
    S_va = state[_IDX_S_VA]
    S_bu = state[_IDX_S_BU]
    S_pro = state[_IDX_S_PRO]
    S_ac = state[_IDX_S_AC]
    S_h2 = state[_IDX_S_H2]
    S_ch4 = state[_IDX_S_CH4]
    S_co2 = state[_IDX_S_CO2]
    S_nh4 = state[_IDX_S_NH4]
    S_I = state[_IDX_S_I]

    X_PS_ch = state[_IDX_X_PS_CH]
    X_PS_pr = state[_IDX_X_PS_PR]
    X_PS_li = state[_IDX_X_PS_LI]
    X_PF_ch = state[_IDX_X_PF_CH]
    X_PF_pr = state[_IDX_X_PF_PR]
    X_PF_li = state[_IDX_X_PF_LI]
    X_S_ch = state[_IDX_X_S_CH]
    X_S_pr = state[_IDX_X_S_PR]
    X_S_li = state[_IDX_X_S_LI]
    X_I = state[_IDX_X_I]

    X_su = state[_IDX_X_SU]
    X_aa = state[_IDX_X_AA]
    X_fa = state[_IDX_X_FA]
    X_c4 = state[_IDX_X_C4]
    X_pro = state[_IDX_X_PRO]
    X_ac = state[_IDX_X_AC]
    X_h2 = state[_IDX_X_H2]

    S_cation = state[_IDX_S_CATION]
    S_anion = state[_IDX_S_ANION]
    S_va_ion = state[_IDX_S_VA_ION]
    S_bu_ion = state[_IDX_S_BU_ION]
    S_pro_ion = state[_IDX_S_PRO_ION]
    S_ac_ion = state[_IDX_S_AC_ION]
    S_hco3 = state[_IDX_S_HCO3]
    S_nh3 = state[_IDX_S_NH3]

    p_gas_h2 = state[_IDX_P_H2]
    p_gas_ch4 = state[_IDX_P_CH4]
    p_gas_co2 = state[_IDX_P_CO2]
    pTOTAL = state[_IDX_P_TOTAL]

    # ---- Influent and hydraulic flow ---------------------------------
    q_ad = float(np.sum(self._Q)) if self._Q is not None else 0.0
    s_in = self._state_input if self._state_input is not None else [0.0] * 37

    # ---- Parameters --------------------------------------------------
    k = self._kinetic
    st = self._stoich
    fr = self._fractions
    ip = self._inhib_params

    f_fa_li = st["f_fa_li"]
    f_ch_bac = st["f_ch_bac"]
    f_pr_bac = st["f_pr_bac"]
    f_li_bac = st["f_li_bac"]
    f_p_bac = st["f_p_bac"]
    fXI_PS = st["fXI_PS"]
    fXI_PF = st["fXI_PF"]
    fSI = st["fSI_hyd"]

    Y_su = k["Y_su"]
    Y_aa = k["Y_aa"]
    Y_fa = k["Y_fa"]
    Y_c4 = k["Y_c4"]
    Y_pro = k["Y_pro"]
    Y_ac = k["Y_ac"]
    Y_h2 = k["Y_h2"]

    # ---- pH (Newton–Raphson charge balance) --------------------------
    S_H = self._calc_ph(
        S_nh4,
        S_nh3,
        S_hco3,
        S_ac_ion,
        S_pro_ion,
        S_bu_ion,
        S_va_ion,
        S_cation,
        S_anion,
        ip["K_w"],
    )
    S_H = max(S_H, 1.0e-14)

    # ---- Inhibition factors ------------------------------------------
    S_IN = S_nh4 + S_nh3
    I_IN = S_IN / (ip["K_S_IN"] + S_IN + 1.0e-20)

    I_pH_aa = self._pH_inhib(S_H, ip["K_pH_aa"], n=1)
    I_pH_fa = self._pH_inhib(S_H, ip["K_pH_aa"], n=2)
    I_pH_c4 = self._pH_inhib(S_H, ip["K_pH_aa"], n=2)
    I_pH_pro = self._pH_inhib(S_H, ip["K_pH_aa"], n=2)
    I_pH_ac = self._pH_inhib(S_H, ip["K_pH_ac"], n=3)
    I_pH_h2 = self._pH_inhib(S_H, ip["K_pH_h2"], n=3)

    I_h2_fa = ip["K_I_h2_fa"] / (ip["K_I_h2_fa"] + S_h2)
    I_h2_c4 = ip["K_I_h2_c4"] / (ip["K_I_h2_c4"] + S_h2)
    I_h2_pro = ip["K_I_h2_pro"] / (ip["K_I_h2_pro"] + S_h2)

    K_nh3_ac_sq = ip["K_I_nh3"] ** 2
    K_nh3_pro_sq = ip["K_I_nh3_pro"] ** 2
    S_nh3_sq = S_nh3 * S_nh3
    I_nh3 = K_nh3_ac_sq / (K_nh3_ac_sq + S_nh3_sq)
    I_nh3_pro = K_nh3_pro_sq / (K_nh3_pro_sq + S_nh3_sq)

    K_co2_h2_sq = ip["K_S_co2_h2"] * ip["K_S_co2_h2"]
    S_co2_sq = S_co2 * S_co2
    I_co2_h2 = S_co2_sq / (K_co2_h2_sq + S_co2_sq + 1.0e-30)

    Ka_pro = ip["K_a_pro"]
    S_HPr = (S_pro / 112.0) * S_H / (S_H + Ka_pro + 1.0e-20)
    I_HPr = ip["K_IH_pro"] / (ip["K_IH_pro"] + S_HPr + 1.0e-20)

    Ka_ac = ip["K_a_ac"]
    S_HAc = (S_ac / 64.0) * S_H / (S_H + Ka_ac + 1.0e-20)
    I_HAc = ip["K_IH_ac"] / (ip["K_IH_ac"] + S_HAc + 1.0e-20)

    I_su = I_pH_aa * I_IN
    I_aa = I_pH_aa * I_IN
    I_fa = I_pH_fa * I_IN * I_h2_fa
    I_c4 = I_pH_c4 * I_IN * I_h2_c4
    I_pro = I_pH_pro * I_IN * I_h2_pro * I_HPr * I_nh3_pro
    I_ac = I_pH_ac * I_IN * I_nh3 * I_HAc
    I_h2 = I_pH_h2 * I_IN * I_co2_h2

    # ---- Process rates -----------------------------------------------
    Rho_dis_PS_ch = k["k_dis_PS"] * X_PS_ch
    Rho_dis_PS_pr = k["k_dis_PS"] * X_PS_pr
    Rho_dis_PS_li = k["k_dis_PS"] * X_PS_li
    Rho_dis_PF_ch = k["k_dis_PF"] * X_PF_ch
    Rho_dis_PF_pr = k["k_dis_PF"] * X_PF_pr
    Rho_dis_PF_li = k["k_dis_PF"] * X_PF_li

    Rho_hyd_ch = k["k_hyd_ch"] * X_S_ch
    Rho_hyd_pr = k["k_hyd_pr"] * X_S_pr
    Rho_hyd_li = k["k_hyd_li"] * X_S_li

    Rho_su = k["k_m_su"] * S_su / (k["K_S_su"] + S_su + 1.0e-20) * X_su * I_su
    Rho_aa = k["k_m_aa"] * S_aa / (k["K_S_aa"] + S_aa + 1.0e-20) * X_aa * I_aa
    Rho_fa = k["k_m_fa"] * S_fa / (k["K_S_fa"] + S_fa + 1.0e-20) * X_fa * I_fa

    S_vbu = S_va + S_bu + 1.0e-20
    Rho_c4_va = k["k_m_c4"] * S_va / (k["K_S_c4"] + S_va + 1.0e-20) * X_c4 * (S_va / S_vbu) * I_c4
    Rho_c4_bu = k["k_m_c4"] * S_bu / (k["K_S_c4"] + S_bu + 1.0e-20) * X_c4 * (S_bu / S_vbu) * I_c4

    Rho_pro = k["k_m_pro"] * S_pro / (k["K_S_pro"] + S_pro + 1.0e-20) * X_pro * I_pro
    Rho_ac = k["k_m_ac"] * S_ac / (k["K_S_ac"] + S_ac + 1.0e-20) * X_ac * I_ac
    Rho_h2 = k["k_m_h2"] * S_h2 / (k["K_S_h2"] + S_h2 + 1.0e-20) * X_h2 * I_h2

    Rho_dec_su = k["k_dec_su"] * X_su
    Rho_dec_aa = k["k_dec_aa"] * X_aa
    Rho_dec_fa = k["k_dec_fa"] * X_fa
    Rho_dec_c4 = k["k_dec_c4"] * X_c4
    Rho_dec_pro = k["k_dec_pro"] * X_pro
    Rho_dec_ac = k["k_dec_ac"] * X_ac
    Rho_dec_h2 = k["k_dec_h2"] * X_h2
    sum_decay = Rho_dec_su + Rho_dec_aa + Rho_dec_fa + Rho_dec_c4 + Rho_dec_pro + Rho_dec_ac + Rho_dec_h2

    # ---- Acid-base rates --------------------------------------------
    k_AB = ip["k_A_B"]
    Rho_A_va = k_AB * (S_va_ion * S_H - ip["K_a_va"] * (S_va - S_va_ion))
    Rho_A_bu = k_AB * (S_bu_ion * S_H - ip["K_a_bu"] * (S_bu - S_bu_ion))
    Rho_A_pro = k_AB * (S_pro_ion * S_H - ip["K_a_pro"] * (S_pro - S_pro_ion))
    Rho_A_ac = k_AB * (S_ac_ion * S_H - ip["K_a_ac"] * (S_ac - S_ac_ion))
    Rho_A_co2 = k_AB * (S_hco3 * S_H - ip["K_a_co2"] * (S_co2 - S_hco3))
    Rho_A_IN = k_AB * (S_nh3 * S_H - ip["K_a_IN"] * (S_nh4 - S_nh3))

    # ---- Gas transfer rates -----------------------------------------
    k_L_a = self._k_L_a
    if self._calibration_params.get("k_L_a") is not None:
        k_L_a = float(self._calibration_params["k_L_a"])

    Rho_T_h2 = k_L_a * (S_h2 - 16.0 * p_gas_h2 / (self._RT * self._K_H_h2)) * (self.V_liq / self._V_gas)
    Rho_T_ch4 = k_L_a * (S_ch4 - 64.0 * p_gas_ch4 / (self._RT * self._K_H_ch4)) * (self.V_liq / self._V_gas)
    S_co2_free = max(S_co2 - S_hco3, 0.0)
    Rho_T_co2 = k_L_a * (S_co2_free - p_gas_co2 / (self._RT * self._K_H_co2)) * (self.V_liq / self._V_gas)

    k_p = self._k_p
    if self._calibration_params.get("k_p") is not None:
        k_p = float(self._calibration_params["k_p"])
    Rho_T_11 = max(
        k_p * (pTOTAL + self._p_gas_h2o - self._p_ext) * (self.V_liq / self._V_gas),
        0.0,
    )

    # ---- Carbon stoichiometry coefficients for S_co2 balance --------
    C = st
    s_hyd_ch = -C["C_ch"] + (1.0 - fSI) * C["C_su"] + fSI * C["C_I_s"]
    s_hyd_pr = -C["C_pr"] + (1.0 - fSI) * C["C_aa"] + fSI * C["C_I_s"]
    s_hyd_li = -C["C_li"] + (1.0 - fSI) * (f_fa_li * C["C_fa"] + (1.0 - f_fa_li) * C["C_su"]) + fSI * C["C_I_s"]

    s_su = (
        -C["C_su"]
        + (1.0 - Y_su) * (fr["f_bu_su"] * C["C_bu"] + fr["f_pro_su"] * C["C_pro"] + fr["f_ac_su"] * C["C_ac"])
        + Y_su * C["C_bac"]
    )
    s_aa = (
        -C["C_aa"]
        + (1.0 - Y_aa)
        * (fr["f_va_aa"] * C["C_va"] + fr["f_bu_aa"] * C["C_bu"] + fr["f_pro_aa"] * C["C_pro"] + fr["f_ac_aa"] * C["C_ac"])
        + Y_aa * C["C_bac"]
    )
    s_fa = -C["C_fa"] + (1.0 - Y_fa) * 0.7 * C["C_ac"] + Y_fa * C["C_bac"]
    s_c4_va = -C["C_va"] + (1.0 - Y_c4) * (0.54 * C["C_pro"] + 0.31 * C["C_ac"]) + Y_c4 * C["C_bac"]
    s_c4_bu = -C["C_bu"] + (1.0 - Y_c4) * 0.8 * C["C_ac"] + Y_c4 * C["C_bac"]
    s_pro = -C["C_pro"] + (1.0 - Y_pro) * 0.57 * C["C_ac"] + Y_pro * C["C_bac"]
    s_ac = -C["C_ac"] + (1.0 - Y_ac) * C["C_ch4"] + Y_ac * C["C_bac"]
    s_h2 = (1.0 - Y_h2) * C["C_ch4"] + Y_h2 * C["C_bac"]
    s_dec = -C["C_bac"] + f_ch_bac * C["C_ch"] + f_pr_bac * C["C_pr"] + f_li_bac * C["C_li"] + f_p_bac * C["C_I_x"]

    Sigma = (
        s_hyd_ch * Rho_hyd_ch
        + s_hyd_pr * Rho_hyd_pr
        + s_hyd_li * Rho_hyd_li
        + s_su * Rho_su
        + s_aa * Rho_aa
        + s_fa * Rho_fa
        + s_c4_va * Rho_c4_va
        + s_c4_bu * Rho_c4_bu
        + s_pro * Rho_pro
        + s_ac * Rho_ac
        + s_h2 * Rho_h2
        + s_dec * sum_decay
    )

    # ---- Differential equations -------------------------------------
    D_in = q_ad / self.V_liq

    # Sludge volume loss (hydrolysis-rate volume balance, after Schlattmann 2011).
    _q_S_loss = self.V_liq * (
        Rho_hyd_ch * (0.9375 / 1550.0) + Rho_hyd_pr * (0.6125 / 1370.0) + Rho_hyd_li * (0.3474 / 920.0)
    )
    self._q_S_loss_last = float(_q_S_loss)

    # If the caller has supplied an outflow (e.g. via a level controller
    # that owns the volume balance), use it. Otherwise enforce volume
    # conservation by setting Q_out = Q_in - q_S_loss.
    if self._Q_out_override is not None:
        _Q_out = max(float(self._Q_out_override), 0.0)
    else:
        _Q_out = max(q_ad - _q_S_loss, 0.0)
    D_out = _Q_out / self.V_liq

    # --- Dissolved (0–11) ---
    diff_S_su = (
        D_in * s_in[0] - D_out * S_su + (1.0 - fSI) * Rho_hyd_ch + (1.0 - fSI) * (1.0 - f_fa_li) * Rho_hyd_li - Rho_su
    )
    diff_S_aa = D_in * s_in[1] - D_out * S_aa + (1.0 - fSI) * Rho_hyd_pr - Rho_aa
    diff_S_fa = D_in * s_in[2] - D_out * S_fa + (1.0 - fSI) * f_fa_li * Rho_hyd_li - Rho_fa
    diff_S_va = D_in * s_in[3] - D_out * S_va + (1.0 - Y_aa) * fr["f_va_aa"] * Rho_aa - Rho_c4_va
    diff_S_bu = (
        D_in * s_in[4]
        - D_out * S_bu
        + (1.0 - Y_su) * fr["f_bu_su"] * Rho_su
        + (1.0 - Y_aa) * fr["f_bu_aa"] * Rho_aa
        - Rho_c4_bu
    )
    diff_S_pro = (
        D_in * s_in[5]
        - D_out * S_pro
        + (1.0 - Y_su) * fr["f_pro_su"] * Rho_su
        + (1.0 - Y_aa) * fr["f_pro_aa"] * Rho_aa
        + (1.0 - Y_c4) * 0.54 * Rho_c4_va
        - Rho_pro
    )
    diff_S_ac = (
        D_in * s_in[6]
        - D_out * S_ac
        + (1.0 - Y_su) * fr["f_ac_su"] * Rho_su
        + (1.0 - Y_aa) * fr["f_ac_aa"] * Rho_aa
        + (1.0 - Y_fa) * 0.7 * Rho_fa
        + (1.0 - Y_c4) * 0.31 * Rho_c4_va
        + (1.0 - Y_c4) * 0.8 * Rho_c4_bu
        + (1.0 - Y_pro) * 0.57 * Rho_pro
        - Rho_ac
    )
    diff_S_h2 = (
        D_in * s_in[7]
        - D_out * S_h2
        + (1.0 - Y_su) * fr["f_h2_su"] * Rho_su
        + (1.0 - Y_aa) * fr["f_h2_aa"] * Rho_aa
        + (1.0 - Y_fa) * 0.3 * Rho_fa
        + (1.0 - Y_c4) * 0.15 * Rho_c4_va
        + (1.0 - Y_c4) * 0.2 * Rho_c4_bu
        + (1.0 - Y_pro) * 0.43 * Rho_pro
        - Rho_h2
        - self._V_gas / self.V_liq * Rho_T_h2
    )
    diff_S_ch4 = (
        D_in * s_in[8]
        - D_out * S_ch4
        + (1.0 - Y_ac) * Rho_ac
        + (1.0 - Y_h2) * Rho_h2
        - self._V_gas / self.V_liq * Rho_T_ch4
    )
    diff_S_co2 = D_in * s_in[9] - D_out * S_co2 - Sigma - self._V_gas / self.V_liq * Rho_T_co2 + Rho_A_co2

    N_bac = st["N_bac"]
    N_aa = st["N_aa"]
    N_I = st["N_I"]
    diff_S_nh4 = (
        D_in * s_in[10]
        - D_out * S_nh4
        - Y_su * N_bac * Rho_su
        + (N_aa - Y_aa * N_bac) * Rho_aa
        - Y_fa * N_bac * Rho_fa
        - Y_c4 * N_bac * Rho_c4_va
        - Y_c4 * N_bac * Rho_c4_bu
        - Y_pro * N_bac * Rho_pro
        - Y_ac * N_bac * Rho_ac
        - Y_h2 * N_bac * Rho_h2
        + (N_bac - f_pr_bac * N_aa - f_p_bac * N_I) * sum_decay
        + Rho_A_IN
    )
    diff_S_I = D_in * s_in[11] - D_out * S_I + fSI * (Rho_hyd_ch + Rho_hyd_pr + Rho_hyd_li)

    # --- Particulate sub-fractions (12–21) ---
    diff_X_PS_ch = D_in * s_in[12] - D_out * X_PS_ch - Rho_dis_PS_ch
    diff_X_PS_pr = D_in * s_in[13] - D_out * X_PS_pr - Rho_dis_PS_pr
    diff_X_PS_li = D_in * s_in[14] - D_out * X_PS_li - Rho_dis_PS_li
    diff_X_PF_ch = D_in * s_in[15] - D_out * X_PF_ch - Rho_dis_PF_ch
    diff_X_PF_pr = D_in * s_in[16] - D_out * X_PF_pr - Rho_dis_PF_pr
    diff_X_PF_li = D_in * s_in[17] - D_out * X_PF_li - Rho_dis_PF_li
    diff_X_S_ch = (
        D_in * s_in[18]
        - D_out * X_S_ch
        + (1.0 - fXI_PS) * Rho_dis_PS_ch
        + (1.0 - fXI_PF) * Rho_dis_PF_ch
        + f_ch_bac * sum_decay
        - Rho_hyd_ch
    )
    diff_X_S_pr = (
        D_in * s_in[19]
        - D_out * X_S_pr
        + (1.0 - fXI_PS) * Rho_dis_PS_pr
        + (1.0 - fXI_PF) * Rho_dis_PF_pr
        + f_pr_bac * sum_decay
        - Rho_hyd_pr
    )
    diff_X_S_li = (
        D_in * s_in[20]
        - D_out * X_S_li
        + (1.0 - fXI_PS) * Rho_dis_PS_li
        + (1.0 - fXI_PF) * Rho_dis_PF_li
        + f_li_bac * sum_decay
        - Rho_hyd_li
    )
    diff_X_I = (
        D_in * s_in[21]
        - D_out * X_I
        + fXI_PS * (Rho_dis_PS_ch + Rho_dis_PS_pr + Rho_dis_PS_li)
        + fXI_PF * (Rho_dis_PF_ch + Rho_dis_PF_pr + Rho_dis_PF_li)
        + f_p_bac * sum_decay
    )

    # --- Biomass (22–28) ---
    diff_X_su = D_in * s_in[22] - D_out * X_su + Y_su * Rho_su - Rho_dec_su
    diff_X_aa = D_in * s_in[23] - D_out * X_aa + Y_aa * Rho_aa - Rho_dec_aa
    diff_X_fa = D_in * s_in[24] - D_out * X_fa + Y_fa * Rho_fa - Rho_dec_fa
    diff_X_c4 = D_in * s_in[25] - D_out * X_c4 + Y_c4 * Rho_c4_va + Y_c4 * Rho_c4_bu - Rho_dec_c4
    diff_X_pro = D_in * s_in[26] - D_out * X_pro + Y_pro * Rho_pro - Rho_dec_pro
    diff_X_ac = D_in * s_in[27] - D_out * X_ac + Y_ac * Rho_ac - Rho_dec_ac
    diff_X_h2 = D_in * s_in[28] - D_out * X_h2 + Y_h2 * Rho_h2 - Rho_dec_h2

    # --- Charge balance (29–36) ---
    diff_S_cation = D_in * s_in[29] - D_out * S_cation
    diff_S_anion = D_in * s_in[30] - D_out * S_anion
    diff_S_va_ion = D_in * s_in[31] - D_out * S_va_ion - Rho_A_va
    diff_S_bu_ion = D_in * s_in[32] - D_out * S_bu_ion - Rho_A_bu
    diff_S_pro_ion = D_in * s_in[33] - D_out * S_pro_ion - Rho_A_pro
    diff_S_ac_ion = D_in * s_in[34] - D_out * S_ac_ion - Rho_A_ac
    diff_S_hco3 = D_in * s_in[35] - D_out * S_hco3 - Rho_A_co2
    diff_S_nh3 = D_in * s_in[36] - D_out * S_nh3 - Rho_A_IN

    # --- Gas phase (37–40) ---
    diff_p_h2 = Rho_T_h2 * self._RT / 16.0 - p_gas_h2 / pTOTAL * Rho_T_11
    diff_p_ch4 = Rho_T_ch4 * self._RT / 64.0 - p_gas_ch4 / pTOTAL * Rho_T_11
    diff_p_co2 = Rho_T_co2 * self._RT - p_gas_co2 / pTOTAL * Rho_T_11
    diff_pTOT = self._RT / 16.0 * Rho_T_h2 + self._RT / 64.0 * Rho_T_ch4 + self._RT * Rho_T_co2 - Rho_T_11

    return (
        diff_S_su,
        diff_S_aa,
        diff_S_fa,
        diff_S_va,
        diff_S_bu,
        diff_S_pro,
        diff_S_ac,
        diff_S_h2,
        diff_S_ch4,
        diff_S_co2,
        diff_S_nh4,
        diff_S_I,
        diff_X_PS_ch,
        diff_X_PS_pr,
        diff_X_PS_li,
        diff_X_PF_ch,
        diff_X_PF_pr,
        diff_X_PF_li,
        diff_X_S_ch,
        diff_X_S_pr,
        diff_X_S_li,
        diff_X_I,
        diff_X_su,
        diff_X_aa,
        diff_X_fa,
        diff_X_c4,
        diff_X_pro,
        diff_X_ac,
        diff_X_h2,
        diff_S_cation,
        diff_S_anion,
        diff_S_va_ion,
        diff_S_bu_ion,
        diff_S_pro_ion,
        diff_S_ac_ion,
        diff_S_hco3,
        diff_S_nh3,
        diff_p_h2,
        diff_p_ch4,
        diff_p_co2,
        diff_pTOT,
    )

__init__(feedstock, V_liq=1977.0, V_gas=304.0, T_ad=308.15)

Initialize the ADM1 model.

Parameters

feedstock : Feedstock Feedstock object (used by create_influent when no external influent DataFrame has been provided via set_influent_dataframe). V_liq : float Liquid digester volume [m³]. V_gas : float Gas headspace volume [m³]. T_ad : float Operating temperature [K] (default 308.15 K = 35 °C).

Source code in pyadm1/core/adm1.py

def __init__(
    self,
    feedstock,
    V_liq: float = 1977.0,
    V_gas: float = 304.0,
    T_ad: float = 308.15,
) -> None:
    """
    Initialize the ADM1 model.

    Parameters
    ----------
    feedstock : Feedstock
        Feedstock object (used by ``create_influent`` when no external
        influent DataFrame has been provided via ``set_influent_dataframe``).
    V_liq : float
        Liquid digester volume [m³].
    V_gas : float
        Gas headspace volume [m³].
    T_ad : float
        Operating temperature [K] (default 308.15 K = 35 °C).
    """
    # --- Reactor volumes ---
    self.V_liq = V_liq
    self._V_gas = V_gas
    self._V_ad = V_liq + V_gas

    # --- Temperature ---
    self._T_ad = T_ad

    # --- Physical constants ---
    self._R = 0.08314  # bar·m³·kmol⁻¹·K⁻¹
    self._T_base = 298.15  # 25 °C reference
    self._p_atm = 1.013

    self._RT = self._R * self._T_ad
    self._p_ext = self._p_atm - 0.0084147 * np.exp(0.054 * (self._T_ad - 273.15))

    # --- Feedstock / influent ---
    self._feedstock = feedstock
    self._Q: Optional[List[float]] = None
    self._state_input: Optional[List[float]] = None

    # --- Calibration overrides ---
    self._calibration_params: dict = {}

    # --- Result-tracking lists ---
    self._Q_GAS: List[float] = []
    self._Q_CH4: List[float] = []
    self._Q_CO2: List[float] = []
    self._Q_H2O: List[float] = []
    self._P_GAS: List[float] = []
    self._pH_l: List[float] = []
    self._FOSTAC: List[float] = []
    self._AcvsPro: List[float] = []
    self._VFA: List[float] = []
    self._TAC: List[float] = []

    # --- Temperature-corrected kinetics (reused across ODE calls) ---
    base_kinetic = ADMParams.get_kinetic_params()
    theta = ADMParams.get_temperature_factors()
    self._kinetic = ADMParams.apply_temperature_corrections(base_kinetic, theta, T_ad)
    # Snapshot of the temperature-corrected defaults so calibration
    # overrides can be reverted cleanly by clear_calibration_parameters().
    self._kinetic_default = dict(self._kinetic)

    # --- Stoichiometric / fraction / inhibition parameters ---
    self._stoich = ADMParams.get_stoichiometric_params()
    self._fractions = ADMParams.get_product_fractions()
    self._inhib_params = ADMParams.get_inhibition_params(self._R, self._T_base, T_ad)

    # --- Gas parameters ---
    p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2 = ADMParams.getADMgasparams(self._R, self._T_base, T_ad)
    self._p_gas_h2o = p_gas_h2o
    self._k_p = k_p
    self._k_L_a = k_L_a
    self._K_H_co2 = K_H_co2
    self._K_H_ch4 = K_H_ch4
    self._K_H_h2 = K_H_h2
    self._K_H_co2_default = K_H_co2
    self._K_H_ch4_default = K_H_ch4
    self._K_H_h2_default = K_H_h2

    # Gas volume reference conditions (theta = 20 °C, 1 atm) — ADM1da convention
    self._T_gas_norm = 293.15
    self._P_gas_norm = 1.01325
    self._NQ = 1000.0 * self._P_gas_norm / (self._R * self._T_gas_norm)

    # Pressure threshold for linearised gas-outlet pipe formula [bar].
    self._dP_min_pipe = 1.0e-3

    # Optional external influent DataFrame.
    self._influent_df: Optional[pd.DataFrame] = None

    # Influent / sludge densities (used by ``create_influent`` only).
    self._rho_in: float = 1000.0
    self._rho_sludge: float = 1000.0

    # When ``_Q_out_override`` is not None, ADM_ODE uses it for the sludge
    # outflow instead of the default ``q_in - q_S_loss`` formula. The
    # caller is then responsible for advancing V_liq between ODE steps.
    # ``_q_S_loss_last`` caches the most recent q_S_loss so the caller
    # can integrate the volume balance externally.
    self._Q_out_override: Optional[float] = None
    self._q_S_loss_last: float = 0.0

calc_gas(pi_Sh2, pi_Sch4, pi_Sco2, pTOTAL)

Calculate biogas production rates, normalised to standard reference conditions (theta = 20 °C, P = 1.01325 bar; ADM1da convention).

Returns

q_gas, q_ch4, q_co2, q_h2o, p_gas

Source code in pyadm1/core/adm1.py

def calc_gas(
    self,
    pi_Sh2: float,
    pi_Sch4: float,
    pi_Sco2: float,
    pTOTAL: float,
) -> Tuple[float, float, float, float, float]:
    """
    Calculate biogas production rates, normalised to standard reference
    conditions (theta = 20 °C, P = 1.01325 bar; ADM1da convention).

    Returns
    -------
    q_gas, q_ch4, q_co2, q_h2o, p_gas
    """
    k_p = self._k_p
    if self._calibration_params.get("k_p") is not None:
        k_p = float(self._calibration_params["k_p"])

    p_total_wet = pTOTAL + self._p_gas_h2o
    q_gas = max(
        k_p * (p_total_wet - self._p_ext) / (self._RT / 1000.0 * self._NQ) * self.V_liq,
        0.0,
    )

    p_gas = pi_Sh2 + pi_Sch4 + pi_Sco2
    p_gas_wet = p_gas + self._p_gas_h2o
    if p_gas_wet > 0.0:
        q_ch4 = max(q_gas * (pi_Sch4 / p_gas_wet), 0.0)
        q_co2 = max(q_gas * (pi_Sco2 / p_gas_wet), 0.0)
        q_h2o = max(q_gas * (self._p_gas_h2o / p_gas_wet), 0.0)
    else:
        q_ch4 = 0.0
        q_co2 = 0.0
        q_h2o = 0.0

    return q_gas, q_ch4, q_co2, q_h2o, p_gas

clear_calibration_parameters()

Clear all calibration parameters and restore defaults.

Source code in pyadm1/core/adm1.py

def clear_calibration_parameters(self) -> None:
    """Clear all calibration parameters and restore defaults."""
    self._calibration_params = {}
    self._kinetic = dict(self._kinetic_default)
    self._K_H_co2 = self._K_H_co2_default
    self._K_H_ch4 = self._K_H_ch4_default
    self._K_H_h2 = self._K_H_h2_default

create_influent(Q, i, rho=None)

Build the influent vector for time step i.

If an external influent DataFrame has been set via set_influent_dataframe(), that DataFrame is used. Otherwise the feedstock object is used to derive the influent.

Source code in pyadm1/core/adm1.py

def create_influent(
    self,
    Q: List[float],
    i: int,
    rho: Optional[List[float]] = None,
) -> None:
    """
    Build the influent vector for time step *i*.

    If an external influent DataFrame has been set via
    ``set_influent_dataframe()``, that DataFrame is used.  Otherwise
    the feedstock object is used to derive the influent.
    """
    if hasattr(self._feedstock, "actual_Q"):
        Q_actual = self._feedstock.actual_Q(Q)
    else:
        Q_actual = list(Q)
    self._Q = Q_actual
    if rho is not None and len(rho) == len(Q_actual):
        q_total = sum(Q_actual)
        if q_total > 0.0:
            self._rho_in = sum(q * r for q, r in zip(Q_actual, rho)) / q_total

    if self._influent_df is not None:
        max_i = len(self._influent_df) - 1
        row = self._influent_df.iloc[min(i, max_i)]
        self._state_input = row[INFLUENT_COLUMNS[:-1]].tolist()
    else:
        df = self._feedstock.get_influent_dataframe(Q)
        max_i = len(df) - 1
        row = df.iloc[min(i, max_i)]
        self._state_input = row[INFLUENT_COLUMNS[:-1]].tolist()

get_calibration_parameters()

Return a copy of the current calibration parameters.

Source code in pyadm1/core/adm1.py

def get_calibration_parameters(self) -> dict:
    """Return a copy of the current calibration parameters."""
    return self._calibration_params.copy()

get_state_size()

Return the number of state variables: 41.

Source code in pyadm1/core/adm1.py

def get_state_size(self) -> int:
    """Return the number of state variables: 41."""
    return STATE_SIZE

print_params_at_current_state(state)

Compute and store process indicators (pH, gas) from the current state.

Source code in pyadm1/core/adm1.py

def print_params_at_current_state(self, state: List[float]) -> None:
    """Compute and store process indicators (pH, gas) from the current state."""
    ip = self._inhib_params
    S_H = self._calc_ph(
        state[_IDX_S_NH4],
        state[_IDX_S_NH3],
        state[_IDX_S_HCO3],
        state[_IDX_S_AC_ION],
        state[_IDX_S_PRO_ION],
        state[_IDX_S_BU_ION],
        state[_IDX_S_VA_ION],
        state[_IDX_S_CATION],
        state[_IDX_S_ANION],
        ip["K_w"],
    )
    pH = -np.log10(max(S_H, 1.0e-14))
    self._track_pH(round(pH, 1))

    q_gas, q_ch4, q_co2, q_h2o, p_gas = self.calc_gas(
        state[_IDX_P_H2],
        state[_IDX_P_CH4],
        state[_IDX_P_CO2],
        state[_IDX_P_TOTAL],
    )
    self._track_gas(q_gas, q_ch4, q_co2, q_h2o, p_gas)

resume_from_broken_simulation(Q_CH4)

Re-populate the methane tracking list after a simulation restart.

Source code in pyadm1/core/adm1.py

def resume_from_broken_simulation(self, Q_CH4: List[float]) -> None:
    """Re-populate the methane tracking list after a simulation restart."""
    for q in Q_CH4:
        self._Q_CH4.append(q)

set_calibration_parameters(parameters)

Set calibration overrides for kinetic / gas parameters.

Kinetic-rate keys (k_m_*, k_hyd_*, k_dis_*, Y_*, K_S_*, k_dec_*) are applied directly into the live kinetic dict consumed by the ODE. Henry-constant keys (K_H_co2, K_H_ch4, K_H_h2) overwrite the corresponding attributes. k_p and k_L_a continue to be looked up from self._calibration_params at their respective call sites.

Source code in pyadm1/core/adm1.py

def set_calibration_parameters(self, parameters: dict) -> None:
    """Set calibration overrides for kinetic / gas parameters.

    Kinetic-rate keys (``k_m_*``, ``k_hyd_*``, ``k_dis_*``, ``Y_*``,
    ``K_S_*``, ``k_dec_*``) are applied directly into the live kinetic
    dict consumed by the ODE. Henry-constant keys (``K_H_co2``,
    ``K_H_ch4``, ``K_H_h2``) overwrite the corresponding attributes.
    ``k_p`` and ``k_L_a`` continue to be looked up from
    ``self._calibration_params`` at their respective call sites.
    """
    self._calibration_params.update(parameters)
    for key, value in parameters.items():
        if key in self._kinetic_default:
            self._kinetic[key] = float(value)
        elif key == "K_H_co2":
            self._K_H_co2 = float(value)
        elif key == "K_H_ch4":
            self._K_H_ch4 = float(value)
        elif key == "K_H_h2":
            self._K_H_h2 = float(value)

set_influent_dataframe(df)

Store an external ADM1 influent DataFrame.

The DataFrame must have columns matching INFLUENT_COLUMNS (37 liquid-state columns + Q). Once set, create_influent() reads from this DataFrame instead of deriving values from the feedstock.

Source code in pyadm1/core/adm1.py

def set_influent_dataframe(self, df: pd.DataFrame) -> None:
    """
    Store an external ADM1 influent DataFrame.

    The DataFrame must have columns matching ``INFLUENT_COLUMNS`` (37
    liquid-state columns + Q).  Once set, ``create_influent()`` reads
    from this DataFrame instead of deriving values from the feedstock.
    """
    missing = [c for c in INFLUENT_COLUMNS if c not in df.columns]
    if missing:
        raise ValueError(f"Influent DataFrame missing columns: {missing}")
    self._influent_df = df.reset_index(drop=True)

set_influent_density(rho_in, rho_sludge=1000.0)

Retained for API compatibility; stored but no longer affects Q_out.

Source code in pyadm1/core/adm1.py

def set_influent_density(self, rho_in: float, rho_sludge: float = 1000.0) -> None:
    """Retained for API compatibility; stored but no longer affects Q_out."""
    self._rho_in = float(rho_in)
    self._rho_sludge = float(rho_sludge)

pyadm1.core.adm1.STATE_SIZE = 41 module-attribute

pyadm1.core.adm1.INFLUENT_COLUMNS = ['S_su', 'S_aa', 'S_fa', 'S_va', 'S_bu', 'S_pro', 'S_ac', 'S_h2', 'S_ch4', 'S_co2', 'S_nh4', 'S_I', 'X_PS_ch', 'X_PS_pr', 'X_PS_li', 'X_PF_ch', 'X_PF_pr', 'X_PF_li', 'X_S_ch', 'X_S_pr', 'X_S_li', 'X_I', 'X_su', 'X_aa', 'X_fa', 'X_c4', 'X_pro', 'X_ac', 'X_h2', 'S_cation', 'S_anion', 'S_va_ion', 'S_bu_ion', 'S_pro_ion', 'S_ac_ion', 'S_hco3_ion', 'S_nh3', 'Q'] module-attribute

pyadm1.core.adm1.get_state_zero_from_csv(csv_file)

Load an ADM1 initial state vector from a CSV file.

The CSV must have columns matching INFLUENT_COLUMNS (without Q) plus the four gas-phase columns (p_gas_h2, p_gas_ch4, p_gas_co2, pTOTAL).

Source code in pyadm1/core/adm1.py

def get_state_zero_from_csv(csv_file: str) -> List[float]:
    """
    Load an ADM1 initial state vector from a CSV file.

    The CSV must have columns matching ``INFLUENT_COLUMNS`` (without ``Q``)
    plus the four gas-phase columns (p_gas_h2, p_gas_ch4, p_gas_co2, pTOTAL).
    """
    df = pd.read_csv(csv_file)
    state = []
    liquid_cols = INFLUENT_COLUMNS[:-1]  # drop Q
    for col in liquid_cols:
        state.append(float(df[col].iloc[0]))
    for col in ("p_gas_h2", "p_gas_ch4", "p_gas_co2", "pTOTAL"):
        state.append(float(df[col].iloc[0]))
    return state

pyadm1.core.adm_params.ADMParams

Static parameter class for the ADM1da model.

Source code in pyadm1/core/adm_params.py

class ADMParams:
    """Static parameter class for the ADM1da model."""

    # *** PUBLIC STATIC GET methods ***

    @staticmethod
    def get_stoichiometric_params() -> dict:
        """
        Return carbon/nitrogen content and disintegration/hydrolysis fractions.

        Returns
        -------
        dict
            Stoichiometric constants (ADM1da defaults; Schlattmann 2011).
        """
        return {
            # --- Carbon content [kmol C / kg COD] ---
            "C_su": 0.0313,
            "C_aa": 0.03,
            "C_fa": 0.0217,
            "C_va": 0.024,
            "C_bu": 0.025,
            "C_pro": 0.0268,
            "C_ac": 0.0313,
            "C_ch4": 0.0156,
            "C_bac": 0.030381,
            "C_ch": 0.0313,
            "C_pr": 0.0306,
            "C_li": 0.022,
            "C_I_s": 0.03,
            "C_I_x": 0.03,
            # --- Nitrogen content [kmol N / kg COD] ---
            "N_bac": 0.005353,
            "N_aa": 0.0076,
            "N_I": 0.06 / 14,
            # --- Lipid hydrolysis fraction ---
            "f_fa_li": 0.95,
            # --- Biomass decay product fractions ---
            "f_ch_bac": 0.2456 * 0.80,
            "f_pr_bac": 0.7093 * 0.80,
            "f_li_bac": 0.0455 * 0.80,
            "f_p_bac": 0.20,
            # --- Disintegration inert fractions ---
            "fXI_PS": 0.0,
            "fXI_PF": 0.0,
            # --- Hydrolysis soluble inert fraction ---
            "fSI_hyd": 0.0,
        }

    @staticmethod
    def get_kinetic_params() -> dict:
        """
        Return kinetic parameters at the reference temperature (35 °C).
        """
        return {
            # --- Disintegration rate constants [d⁻¹] ---
            "k_dis_PS": 0.04,
            "k_dis_PF": 0.4,
            # --- Hydrolysis rate constants [d⁻¹] ---
            "k_hyd_ch": 4.0,
            "k_hyd_pr": 4.0,
            "k_hyd_li": 4.0,
            # --- Maximum uptake rates [d⁻¹] ---
            "k_m_su": 30.0,
            "k_m_aa": 50.0,
            "k_m_fa": 6.0,
            "k_m_c4": 20.0,
            "k_m_pro": 13.0,
            "k_m_ac": 8.0,
            "k_m_h2": 35.0,
            # --- Half-saturation constants [kg COD m⁻³] ---
            "K_S_su": 0.5,
            "K_S_aa": 0.3,
            "K_S_fa": 0.4,
            "K_S_c4": 0.2,
            "K_S_pro": 0.1,
            "K_S_ac": 0.15,
            "K_S_h2": 7.0e-6,
            # --- Per-organism decay rates [d⁻¹] ---
            "k_dec_su": 0.02,
            "k_dec_aa": 0.02,
            "k_dec_fa": 0.02,
            "k_dec_c4": 0.02,
            "k_dec_pro": 0.02,
            "k_dec_ac": 0.04,
            "k_dec_h2": 0.02,
            # --- Yield coefficients ---
            "Y_su": 0.10,
            "Y_aa": 0.08,
            "Y_fa": 0.06,
            "Y_c4": 0.06,
            "Y_pro": 0.04,
            "Y_ac": 0.05,
            "Y_h2": 0.06,
        }

    @staticmethod
    def get_temperature_factors() -> dict:
        """
        Return Arrhenius θ correction factors per organism/process group.
        """
        return {
            "theta_dis": np.exp(0.024),
            "theta_hyd": np.exp(0.024),
            "theta_su": np.exp(0.069),
            "theta_aa": np.exp(0.069),
            "theta_fa": np.exp(0.055),
            "theta_c4": np.exp(0.055),
            "theta_pro": np.exp(0.055),
            "theta_ac": np.exp(0.055),
            "theta_h2": np.exp(0.069),
            "theta_dec_su_aa_h2": np.exp(0.069),
            "theta_dec_fa_c4_pro_ac": np.exp(0.055),
        }

    @staticmethod
    def get_product_fractions() -> dict:
        """Return fermentation product fractions."""
        return {
            "f_h2_su": 0.19,
            "f_bu_su": 0.13,
            "f_pro_su": 0.27,
            "f_ac_su": 0.41,
            "f_h2_aa": 0.06,
            "f_va_aa": 0.23,
            "f_bu_aa": 0.26,
            "f_pro_aa": 0.05,
            "f_ac_aa": 0.40,
        }

    @staticmethod
    def get_inhibition_params(R: float, T_base: float, T_ad: float) -> dict:
        """
        Return inhibition-related parameters and acid-base constants.

        Parameters
        ----------
        R : float
            Gas constant [bar m³ kmol⁻¹ K⁻¹].
        T_base : float
            Reference temperature [K] (25 °C = 298.15 K).
        T_ad : float
            Operating temperature [K].
        """
        pH_LL_aa, pH_UL_aa = 4.0, 5.5
        pH_LL_ac, pH_UL_ac = 6.0, 7.0
        pH_LL_h2, pH_UL_h2 = 5.0, 6.0

        K_pH_aa = 10.0 ** (-(pH_LL_aa + pH_UL_aa) / 2.0)
        K_pH_ac = 10.0 ** (-(pH_LL_ac + pH_UL_ac) / 2.0)
        K_pH_h2 = 10.0 ** (-(pH_LL_h2 + pH_UL_h2) / 2.0)

        K_w = 10.0**-14.0 * np.exp((55900.0 / (100.0 * R)) * (1.0 / T_base - 1.0 / T_ad))
        K_a_va = 10.0**-4.86
        K_a_bu = 10.0**-4.82
        K_a_pro = 10.0**-4.88
        K_a_ac = 10.0**-4.76
        K_a_co2 = 10.0**-6.35 * np.exp((7646.0 / (100.0 * R)) * (1.0 / T_base - 1.0 / T_ad))
        K_a_IN = 10.0**-9.25 * np.exp((51965.0 / (100.0 * R)) * (1.0 / T_base - 1.0 / T_ad))

        dT_C = T_ad - 308.15

        return {
            "K_pH_aa": K_pH_aa,
            "K_pH_ac": K_pH_ac,
            "K_pH_h2": K_pH_h2,
            "K_w": K_w,
            "K_a_va": K_a_va,
            "K_a_bu": K_a_bu,
            "K_a_pro": K_a_pro,
            "K_a_ac": K_a_ac,
            "K_a_co2": K_a_co2,
            "K_a_IN": K_a_IN,
            "k_A_B": 1.0e8,
            "K_S_IN": 1.0e-4,
            "K_I_nh3": 0.0018 * np.exp(0.086 * dT_C),
            "K_I_h2_fa": 5.0e-6 * np.exp(0.080 * dT_C),
            "K_I_h2_c4": 1.0e-5 * np.exp(0.080 * dT_C),
            "K_I_h2_pro": 3.5e-6 * np.exp(0.080 * dT_C),
            "K_IH_pro": 8.0e-4,
            "K_IH_ac": 2.417e-3,
            "K_I_nh3_pro": 0.0019 * np.exp(0.060 * dT_C),
            "K_S_co2_h2": 5.0e-5,
            "K_I_ac_xfa": 4.0 * np.exp(0.080 * dT_C),
            "K_I_ac_xc4": 4.0 * np.exp(0.080 * dT_C),
        }

    @staticmethod
    def apply_temperature_corrections(kinetic: dict, theta: dict, T_ad: float) -> dict:
        """
        Return a copy of *kinetic* with rates corrected to *T_ad*.

        Uses: k(T) = k(35 °C) · θ^(T[°C] − 35)
        """
        dT = T_ad - 308.15
        corrected = dict(kinetic)

        corrected["k_dis_PS"] = kinetic["k_dis_PS"] * theta["theta_dis"] ** dT
        corrected["k_dis_PF"] = kinetic["k_dis_PF"] * theta["theta_dis"] ** dT
        corrected["k_hyd_ch"] = kinetic["k_hyd_ch"] * theta["theta_hyd"] ** dT
        corrected["k_hyd_pr"] = kinetic["k_hyd_pr"] * theta["theta_hyd"] ** dT
        corrected["k_hyd_li"] = kinetic["k_hyd_li"] * theta["theta_hyd"] ** dT
        corrected["k_m_su"] = kinetic["k_m_su"] * theta["theta_su"] ** dT
        corrected["k_m_aa"] = kinetic["k_m_aa"] * theta["theta_aa"] ** dT
        corrected["k_m_fa"] = kinetic["k_m_fa"] * theta["theta_fa"] ** dT
        corrected["k_m_c4"] = kinetic["k_m_c4"] * theta["theta_c4"] ** dT
        corrected["k_m_pro"] = kinetic["k_m_pro"] * theta["theta_pro"] ** dT
        corrected["k_m_ac"] = kinetic["k_m_ac"] * theta["theta_ac"] ** dT
        corrected["k_m_h2"] = kinetic["k_m_h2"] * theta["theta_h2"] ** dT
        for dec_key in ("k_dec_su", "k_dec_aa", "k_dec_h2"):
            corrected[dec_key] = kinetic[dec_key] * theta["theta_dec_su_aa_h2"] ** dT
        for dec_key in ("k_dec_fa", "k_dec_c4", "k_dec_pro", "k_dec_ac"):
            corrected[dec_key] = kinetic[dec_key] * theta["theta_dec_fa_c4_pro_ac"] ** dT

        return corrected

    @staticmethod
    def getADMgasparams(R: float, T_base: float, T_ad: float) -> Tuple[float, float, float, float, float, float]:
        """
        Get gas phase parameters including Henry constants.

        Parameters
        ----------
        R : float
            Gas constant [bar·m³·kmol⁻¹·K⁻¹]
        T_base : float
            Reference temperature [K] (298.15)
        T_ad : float
            Digester temperature [K]

        Returns
        -------
        Tuple[float, float, float, float, float, float]
            p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2
        """
        p_gas_h2o = 0.0313 * np.exp(5290 * (1 / T_base - 1 / T_ad))
        k_p = 1.0e4
        k_L_a = 200.0

        # Henry's law solubility coefficients with van't Hoff T-correction.
        T_ref_H = 308.15
        R_J = 100.0 * R
        dH_co2, dH_ch4, dH_h2 = 19410.0, 14240.0, 4180.0
        H_co2 = 0.0271 * np.exp(-dH_co2 / R_J * (1.0 / T_ref_H - 1.0 / T_ad))
        H_ch4 = 0.00116 * np.exp(-dH_ch4 / R_J * (1.0 / T_ref_H - 1.0 / T_ad))
        H_h2 = 7.38e-4 * np.exp(-dH_h2 / R_J * (1.0 / T_ref_H - 1.0 / T_ad))
        K_H_co2 = 1 / (H_co2 * R * T_ad)
        K_H_ch4 = 1 / (H_ch4 * R * T_ad)
        K_H_h2 = 1 / (H_h2 * R * T_ad)

        return p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2

Functions

apply_temperature_corrections(kinetic, theta, T_ad) staticmethod

Return a copy of kinetic with rates corrected to T_ad.

Uses: k(T) = k(35 °C) · θ^(T[°C] − 35)

Source code in pyadm1/core/adm_params.py

@staticmethod
def apply_temperature_corrections(kinetic: dict, theta: dict, T_ad: float) -> dict:
    """
    Return a copy of *kinetic* with rates corrected to *T_ad*.

    Uses: k(T) = k(35 °C) · θ^(T[°C] − 35)
    """
    dT = T_ad - 308.15
    corrected = dict(kinetic)

    corrected["k_dis_PS"] = kinetic["k_dis_PS"] * theta["theta_dis"] ** dT
    corrected["k_dis_PF"] = kinetic["k_dis_PF"] * theta["theta_dis"] ** dT
    corrected["k_hyd_ch"] = kinetic["k_hyd_ch"] * theta["theta_hyd"] ** dT
    corrected["k_hyd_pr"] = kinetic["k_hyd_pr"] * theta["theta_hyd"] ** dT
    corrected["k_hyd_li"] = kinetic["k_hyd_li"] * theta["theta_hyd"] ** dT
    corrected["k_m_su"] = kinetic["k_m_su"] * theta["theta_su"] ** dT
    corrected["k_m_aa"] = kinetic["k_m_aa"] * theta["theta_aa"] ** dT
    corrected["k_m_fa"] = kinetic["k_m_fa"] * theta["theta_fa"] ** dT
    corrected["k_m_c4"] = kinetic["k_m_c4"] * theta["theta_c4"] ** dT
    corrected["k_m_pro"] = kinetic["k_m_pro"] * theta["theta_pro"] ** dT
    corrected["k_m_ac"] = kinetic["k_m_ac"] * theta["theta_ac"] ** dT
    corrected["k_m_h2"] = kinetic["k_m_h2"] * theta["theta_h2"] ** dT
    for dec_key in ("k_dec_su", "k_dec_aa", "k_dec_h2"):
        corrected[dec_key] = kinetic[dec_key] * theta["theta_dec_su_aa_h2"] ** dT
    for dec_key in ("k_dec_fa", "k_dec_c4", "k_dec_pro", "k_dec_ac"):
        corrected[dec_key] = kinetic[dec_key] * theta["theta_dec_fa_c4_pro_ac"] ** dT

    return corrected

getADMgasparams(R, T_base, T_ad) staticmethod

Get gas phase parameters including Henry constants.

Parameters

R : float Gas constant [bar·m³·kmol⁻¹·K⁻¹] T_base : float Reference temperature [K] (298.15) T_ad : float Digester temperature [K]

Returns

Tuple[float, float, float, float, float, float] p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2

Source code in pyadm1/core/adm_params.py

@staticmethod
def getADMgasparams(R: float, T_base: float, T_ad: float) -> Tuple[float, float, float, float, float, float]:
    """
    Get gas phase parameters including Henry constants.

    Parameters
    ----------
    R : float
        Gas constant [bar·m³·kmol⁻¹·K⁻¹]
    T_base : float
        Reference temperature [K] (298.15)
    T_ad : float
        Digester temperature [K]

    Returns
    -------
    Tuple[float, float, float, float, float, float]
        p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2
    """
    p_gas_h2o = 0.0313 * np.exp(5290 * (1 / T_base - 1 / T_ad))
    k_p = 1.0e4
    k_L_a = 200.0

    # Henry's law solubility coefficients with van't Hoff T-correction.
    T_ref_H = 308.15
    R_J = 100.0 * R
    dH_co2, dH_ch4, dH_h2 = 19410.0, 14240.0, 4180.0
    H_co2 = 0.0271 * np.exp(-dH_co2 / R_J * (1.0 / T_ref_H - 1.0 / T_ad))
    H_ch4 = 0.00116 * np.exp(-dH_ch4 / R_J * (1.0 / T_ref_H - 1.0 / T_ad))
    H_h2 = 7.38e-4 * np.exp(-dH_h2 / R_J * (1.0 / T_ref_H - 1.0 / T_ad))
    K_H_co2 = 1 / (H_co2 * R * T_ad)
    K_H_ch4 = 1 / (H_ch4 * R * T_ad)
    K_H_h2 = 1 / (H_h2 * R * T_ad)

    return p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2

get_inhibition_params(R, T_base, T_ad) staticmethod

Return inhibition-related parameters and acid-base constants.

Parameters

R : float Gas constant [bar m³ kmol⁻¹ K⁻¹]. T_base : float Reference temperature [K] (25 °C = 298.15 K). T_ad : float Operating temperature [K].

Source code in pyadm1/core/adm_params.py

@staticmethod
def get_inhibition_params(R: float, T_base: float, T_ad: float) -> dict:
    """
    Return inhibition-related parameters and acid-base constants.

    Parameters
    ----------
    R : float
        Gas constant [bar m³ kmol⁻¹ K⁻¹].
    T_base : float
        Reference temperature [K] (25 °C = 298.15 K).
    T_ad : float
        Operating temperature [K].
    """
    pH_LL_aa, pH_UL_aa = 4.0, 5.5
    pH_LL_ac, pH_UL_ac = 6.0, 7.0
    pH_LL_h2, pH_UL_h2 = 5.0, 6.0

    K_pH_aa = 10.0 ** (-(pH_LL_aa + pH_UL_aa) / 2.0)
    K_pH_ac = 10.0 ** (-(pH_LL_ac + pH_UL_ac) / 2.0)
    K_pH_h2 = 10.0 ** (-(pH_LL_h2 + pH_UL_h2) / 2.0)

    K_w = 10.0**-14.0 * np.exp((55900.0 / (100.0 * R)) * (1.0 / T_base - 1.0 / T_ad))
    K_a_va = 10.0**-4.86
    K_a_bu = 10.0**-4.82
    K_a_pro = 10.0**-4.88
    K_a_ac = 10.0**-4.76
    K_a_co2 = 10.0**-6.35 * np.exp((7646.0 / (100.0 * R)) * (1.0 / T_base - 1.0 / T_ad))
    K_a_IN = 10.0**-9.25 * np.exp((51965.0 / (100.0 * R)) * (1.0 / T_base - 1.0 / T_ad))

    dT_C = T_ad - 308.15

    return {
        "K_pH_aa": K_pH_aa,
        "K_pH_ac": K_pH_ac,
        "K_pH_h2": K_pH_h2,
        "K_w": K_w,
        "K_a_va": K_a_va,
        "K_a_bu": K_a_bu,
        "K_a_pro": K_a_pro,
        "K_a_ac": K_a_ac,
        "K_a_co2": K_a_co2,
        "K_a_IN": K_a_IN,
        "k_A_B": 1.0e8,
        "K_S_IN": 1.0e-4,
        "K_I_nh3": 0.0018 * np.exp(0.086 * dT_C),
        "K_I_h2_fa": 5.0e-6 * np.exp(0.080 * dT_C),
        "K_I_h2_c4": 1.0e-5 * np.exp(0.080 * dT_C),
        "K_I_h2_pro": 3.5e-6 * np.exp(0.080 * dT_C),
        "K_IH_pro": 8.0e-4,
        "K_IH_ac": 2.417e-3,
        "K_I_nh3_pro": 0.0019 * np.exp(0.060 * dT_C),
        "K_S_co2_h2": 5.0e-5,
        "K_I_ac_xfa": 4.0 * np.exp(0.080 * dT_C),
        "K_I_ac_xc4": 4.0 * np.exp(0.080 * dT_C),
    }

get_kinetic_params() staticmethod

Return kinetic parameters at the reference temperature (35 °C).

Source code in pyadm1/core/adm_params.py

@staticmethod
def get_kinetic_params() -> dict:
    """
    Return kinetic parameters at the reference temperature (35 °C).
    """
    return {
        # --- Disintegration rate constants [d⁻¹] ---
        "k_dis_PS": 0.04,
        "k_dis_PF": 0.4,
        # --- Hydrolysis rate constants [d⁻¹] ---
        "k_hyd_ch": 4.0,
        "k_hyd_pr": 4.0,
        "k_hyd_li": 4.0,
        # --- Maximum uptake rates [d⁻¹] ---
        "k_m_su": 30.0,
        "k_m_aa": 50.0,
        "k_m_fa": 6.0,
        "k_m_c4": 20.0,
        "k_m_pro": 13.0,
        "k_m_ac": 8.0,
        "k_m_h2": 35.0,
        # --- Half-saturation constants [kg COD m⁻³] ---
        "K_S_su": 0.5,
        "K_S_aa": 0.3,
        "K_S_fa": 0.4,
        "K_S_c4": 0.2,
        "K_S_pro": 0.1,
        "K_S_ac": 0.15,
        "K_S_h2": 7.0e-6,
        # --- Per-organism decay rates [d⁻¹] ---
        "k_dec_su": 0.02,
        "k_dec_aa": 0.02,
        "k_dec_fa": 0.02,
        "k_dec_c4": 0.02,
        "k_dec_pro": 0.02,
        "k_dec_ac": 0.04,
        "k_dec_h2": 0.02,
        # --- Yield coefficients ---
        "Y_su": 0.10,
        "Y_aa": 0.08,
        "Y_fa": 0.06,
        "Y_c4": 0.06,
        "Y_pro": 0.04,
        "Y_ac": 0.05,
        "Y_h2": 0.06,
    }

get_product_fractions() staticmethod

Return fermentation product fractions.

Source code in pyadm1/core/adm_params.py

@staticmethod
def get_product_fractions() -> dict:
    """Return fermentation product fractions."""
    return {
        "f_h2_su": 0.19,
        "f_bu_su": 0.13,
        "f_pro_su": 0.27,
        "f_ac_su": 0.41,
        "f_h2_aa": 0.06,
        "f_va_aa": 0.23,
        "f_bu_aa": 0.26,
        "f_pro_aa": 0.05,
        "f_ac_aa": 0.40,
    }

get_stoichiometric_params() staticmethod

Return carbon/nitrogen content and disintegration/hydrolysis fractions.

Returns

dict Stoichiometric constants (ADM1da defaults; Schlattmann 2011).

Source code in pyadm1/core/adm_params.py

@staticmethod
def get_stoichiometric_params() -> dict:
    """
    Return carbon/nitrogen content and disintegration/hydrolysis fractions.

    Returns
    -------
    dict
        Stoichiometric constants (ADM1da defaults; Schlattmann 2011).
    """
    return {
        # --- Carbon content [kmol C / kg COD] ---
        "C_su": 0.0313,
        "C_aa": 0.03,
        "C_fa": 0.0217,
        "C_va": 0.024,
        "C_bu": 0.025,
        "C_pro": 0.0268,
        "C_ac": 0.0313,
        "C_ch4": 0.0156,
        "C_bac": 0.030381,
        "C_ch": 0.0313,
        "C_pr": 0.0306,
        "C_li": 0.022,
        "C_I_s": 0.03,
        "C_I_x": 0.03,
        # --- Nitrogen content [kmol N / kg COD] ---
        "N_bac": 0.005353,
        "N_aa": 0.0076,
        "N_I": 0.06 / 14,
        # --- Lipid hydrolysis fraction ---
        "f_fa_li": 0.95,
        # --- Biomass decay product fractions ---
        "f_ch_bac": 0.2456 * 0.80,
        "f_pr_bac": 0.7093 * 0.80,
        "f_li_bac": 0.0455 * 0.80,
        "f_p_bac": 0.20,
        # --- Disintegration inert fractions ---
        "fXI_PS": 0.0,
        "fXI_PF": 0.0,
        # --- Hydrolysis soluble inert fraction ---
        "fSI_hyd": 0.0,
    }

get_temperature_factors() staticmethod

Return Arrhenius θ correction factors per organism/process group.

Source code in pyadm1/core/adm_params.py

@staticmethod
def get_temperature_factors() -> dict:
    """
    Return Arrhenius θ correction factors per organism/process group.
    """
    return {
        "theta_dis": np.exp(0.024),
        "theta_hyd": np.exp(0.024),
        "theta_su": np.exp(0.069),
        "theta_aa": np.exp(0.069),
        "theta_fa": np.exp(0.055),
        "theta_c4": np.exp(0.055),
        "theta_pro": np.exp(0.055),
        "theta_ac": np.exp(0.055),
        "theta_h2": np.exp(0.069),
        "theta_dec_su_aa_h2": np.exp(0.069),
        "theta_dec_fa_c4_pro_ac": np.exp(0.055),
    }

pyadm1.core.solver.ODESolver

ODE solver wrapper for ADM1 system.

Provides a clean interface to scipy's solve_ivp with appropriate settings for stiff biogas process ODEs. Uses BDF (Backward Differentiation Formula) method which is suitable for stiff systems.

Example

def ode_func(t, y): ... return [-0.5 * y[0], 0.5 * y[0] - 0.1 * y[1]] solver = ODESolver() result = solver.solve(ode_func, [0, 10], [1.0, 0.0]) print(result.y[:, -1]) # Final state

Source code in pyadm1/core/solver.py

class ODESolver:
    """
    ODE solver wrapper for ADM1 system.

    Provides a clean interface to scipy's solve_ivp with appropriate settings
    for stiff biogas process ODEs. Uses BDF (Backward Differentiation Formula)
    method which is suitable for stiff systems.

    Example:
        >>> def ode_func(t, y):
        ...     return [-0.5 * y[0], 0.5 * y[0] - 0.1 * y[1]]
        >>> solver = ODESolver()
        >>> result = solver.solve(ode_func, [0, 10], [1.0, 0.0])
        >>> print(result.y[:, -1])  # Final state
    """

    def __init__(self, config: Optional[SolverConfig] = None):
        """
        Initialize ODE solver with configuration.

        Args:
            config: Solver configuration. If None, uses default BDF settings.
        """
        self.config = config or SolverConfig()

    def solve(
        self,
        fun: Callable[[float, List[float]], List[float]],
        t_span: Tuple[float, float],
        y0: List[float],
        t_eval: Optional[np.ndarray] = None,
        dense_output: bool = False,
    ) -> "OdeResult":
        """
        Solve ODE system over time span.

        Args:
            fun: Right-hand side of ODE system dy/dt = fun(t, y)
            t_span: Integration time span (t_start, t_end) [days]
            y0: Initial state vector
            t_eval: Times at which to store solution. If None, uses automatic
                time points with 0.05 day resolution
            dense_output: If True, returns continuous solution object

        Returns:
            OdeResult object with solution (from scipy.integrate.solve_ivp)

        Raises:
            RuntimeError: If integration fails
        """
        # Set default evaluation times if not provided
        if t_eval is None:
            t_start, t_end = float(t_span[0]), float(t_span[1])
            step = 0.05
            t_eval = np.arange(t_start, t_end, step, dtype=float)

            # Always include both interval bounds to avoid returning only t_start
            # when (t_end - t_start) < step.
            if t_eval.size == 0 or not np.isclose(t_eval[0], t_start):
                t_eval = np.insert(t_eval, 0, t_start)
            if not np.isclose(t_eval[-1], t_end):
                t_eval = np.append(t_eval, t_end)

        # Prepare solver arguments
        solver_args = {
            "method": self.config.method,
            "rtol": self.config.rtol,
            "atol": self.config.atol,
            "dense_output": dense_output,
        }

        # Add optional step size constraints
        # scipy.solve_ivp only accepts min_step for LSODA.
        if self.config.min_step is not None and self.config.method.upper() == "LSODA":
            solver_args["min_step"] = self.config.min_step
        if self.config.max_step is not None:
            solver_args["max_step"] = self.config.max_step
        if self.config.first_step is not None:
            solver_args["first_step"] = self.config.first_step

        try:
            result = scipy.integrate.solve_ivp(fun=fun, t_span=t_span, y0=y0, t_eval=t_eval, **solver_args)

            if not result.success:
                raise RuntimeError(f"ODE integration failed: {result.message}")

            return result

        except Exception as e:
            raise RuntimeError(f"Error during ODE integration: {str(e)}") from e

    def solve_to_steady_state(
        self,
        fun: Callable[[float, List[float]], List[float]],
        y0: List[float],
        max_time: float = 1000.0,
        steady_state_tol: float = 1e-6,
        check_interval: float = 10.0,
    ) -> Tuple[List[float], float, bool]:
        """
        Integrate until steady state is reached or max time exceeded.

        Args:
            fun: Right-hand side of ODE system
            y0: Initial state vector
            max_time: Maximum integration time [days]
            steady_state_tol: Tolerance for steady state detection
            check_interval: Interval for checking steady state [days]

        Returns:
            Tuple of (final_state, final_time, converged)
            - final_state: State vector at end of integration
            - final_time: Time at end of integration [days]
            - converged: True if steady state was reached
        """
        current_time = 0.0
        current_state = y0

        while current_time < max_time:
            # Integrate for check_interval
            t_span = (current_time, current_time + check_interval)
            result = self.solve(fun, t_span, current_state)

            # Get final state
            new_state = result.y[:, -1]

            # Check if steady state reached
            state_change = np.linalg.norm(new_state - current_state)
            state_norm = np.linalg.norm(new_state)

            if state_norm > 0:
                relative_change = state_change / state_norm
                if relative_change < steady_state_tol:
                    return list(new_state), current_time + check_interval, True

            # Update for next iteration
            current_state = new_state
            current_time += check_interval

        # Max time exceeded without reaching steady state
        return list(current_state), current_time, False

    def solve_sequential(
        self,
        fun: Callable[[float, List[float]], List[float]],
        t_points: List[float],
        y0: List[float],
    ) -> List[List[float]]:
        """
        Solve ODE system sequentially through multiple time points.

        Useful for simulations where conditions change at specific times
        (e.g., substrate feed changes).

        Args:
            fun: Right-hand side of ODE system
            t_points: List of time points [days]
            y0: Initial state vector

        Returns:
            List of state vectors at each time point
        """
        states = [y0]
        current_state = y0

        for i in range(len(t_points) - 1):
            t_span = (t_points[i], t_points[i + 1])
            result = self.solve(fun, t_span, current_state)
            current_state = result.y[:, -1].tolist()
            states.append(current_state)

        return states

Functions

__init__(config=None)

Initialize ODE solver with configuration.

Parameters:

Name	Type	Description	Default
config	Optional[SolverConfig]	Solver configuration. If None, uses default BDF settings.	None

Source code in pyadm1/core/solver.py

def __init__(self, config: Optional[SolverConfig] = None):
    """
    Initialize ODE solver with configuration.

    Args:
        config: Solver configuration. If None, uses default BDF settings.
    """
    self.config = config or SolverConfig()

solve(fun, t_span, y0, t_eval=None, dense_output=False)

Solve ODE system over time span.

Parameters:

Name	Type	Description	Default
fun	Callable[[float, List[float]], List[float]]	Right-hand side of ODE system dy/dt = fun(t, y)	required
t_span	Tuple[float, float]	Integration time span (t_start, t_end) [days]	required
y0	List[float]	Initial state vector	required
t_eval	Optional[ndarray]	Times at which to store solution. If None, uses automatic time points with 0.05 day resolution	None
dense_output	bool	If True, returns continuous solution object	False

Returns:

Type	Description
OdeResult	OdeResult object with solution (from scipy.integrate.solve_ivp)

Raises:

Type	Description
RuntimeError	If integration fails

Source code in pyadm1/core/solver.py

def solve(
    self,
    fun: Callable[[float, List[float]], List[float]],
    t_span: Tuple[float, float],
    y0: List[float],
    t_eval: Optional[np.ndarray] = None,
    dense_output: bool = False,
) -> "OdeResult":
    """
    Solve ODE system over time span.

    Args:
        fun: Right-hand side of ODE system dy/dt = fun(t, y)
        t_span: Integration time span (t_start, t_end) [days]
        y0: Initial state vector
        t_eval: Times at which to store solution. If None, uses automatic
            time points with 0.05 day resolution
        dense_output: If True, returns continuous solution object

    Returns:
        OdeResult object with solution (from scipy.integrate.solve_ivp)

    Raises:
        RuntimeError: If integration fails
    """
    # Set default evaluation times if not provided
    if t_eval is None:
        t_start, t_end = float(t_span[0]), float(t_span[1])
        step = 0.05
        t_eval = np.arange(t_start, t_end, step, dtype=float)

        # Always include both interval bounds to avoid returning only t_start
        # when (t_end - t_start) < step.
        if t_eval.size == 0 or not np.isclose(t_eval[0], t_start):
            t_eval = np.insert(t_eval, 0, t_start)
        if not np.isclose(t_eval[-1], t_end):
            t_eval = np.append(t_eval, t_end)

    # Prepare solver arguments
    solver_args = {
        "method": self.config.method,
        "rtol": self.config.rtol,
        "atol": self.config.atol,
        "dense_output": dense_output,
    }

    # Add optional step size constraints
    # scipy.solve_ivp only accepts min_step for LSODA.
    if self.config.min_step is not None and self.config.method.upper() == "LSODA":
        solver_args["min_step"] = self.config.min_step
    if self.config.max_step is not None:
        solver_args["max_step"] = self.config.max_step
    if self.config.first_step is not None:
        solver_args["first_step"] = self.config.first_step

    try:
        result = scipy.integrate.solve_ivp(fun=fun, t_span=t_span, y0=y0, t_eval=t_eval, **solver_args)

        if not result.success:
            raise RuntimeError(f"ODE integration failed: {result.message}")

        return result

    except Exception as e:
        raise RuntimeError(f"Error during ODE integration: {str(e)}") from e

solve_sequential(fun, t_points, y0)

Solve ODE system sequentially through multiple time points.

Useful for simulations where conditions change at specific times (e.g., substrate feed changes).

Parameters:

Name	Type	Description	Default
fun	Callable[[float, List[float]], List[float]]	Right-hand side of ODE system	required
t_points	List[float]	List of time points [days]	required
y0	List[float]	Initial state vector	required

Returns:

Type	Description
List[List[float]]	List of state vectors at each time point

Source code in pyadm1/core/solver.py

def solve_sequential(
    self,
    fun: Callable[[float, List[float]], List[float]],
    t_points: List[float],
    y0: List[float],
) -> List[List[float]]:
    """
    Solve ODE system sequentially through multiple time points.

    Useful for simulations where conditions change at specific times
    (e.g., substrate feed changes).

    Args:
        fun: Right-hand side of ODE system
        t_points: List of time points [days]
        y0: Initial state vector

    Returns:
        List of state vectors at each time point
    """
    states = [y0]
    current_state = y0

    for i in range(len(t_points) - 1):
        t_span = (t_points[i], t_points[i + 1])
        result = self.solve(fun, t_span, current_state)
        current_state = result.y[:, -1].tolist()
        states.append(current_state)

    return states

solve_to_steady_state(fun, y0, max_time=1000.0, steady_state_tol=1e-06, check_interval=10.0)

Integrate until steady state is reached or max time exceeded.

Parameters:

Name	Type	Description	Default
fun	Callable[[float, List[float]], List[float]]	Right-hand side of ODE system	required
y0	List[float]	Initial state vector	required
max_time	float	Maximum integration time [days]	1000.0
steady_state_tol	float	Tolerance for steady state detection	1e-06
check_interval	float	Interval for checking steady state [days]	10.0

Returns:

Type	Description
List[float]	Tuple of (final_state, final_time, converged)
float	final_state: State vector at end of integration
bool	final_time: Time at end of integration [days]
Tuple[List[float], float, bool]	converged: True if steady state was reached

Source code in pyadm1/core/solver.py

def solve_to_steady_state(
    self,
    fun: Callable[[float, List[float]], List[float]],
    y0: List[float],
    max_time: float = 1000.0,
    steady_state_tol: float = 1e-6,
    check_interval: float = 10.0,
) -> Tuple[List[float], float, bool]:
    """
    Integrate until steady state is reached or max time exceeded.

    Args:
        fun: Right-hand side of ODE system
        y0: Initial state vector
        max_time: Maximum integration time [days]
        steady_state_tol: Tolerance for steady state detection
        check_interval: Interval for checking steady state [days]

    Returns:
        Tuple of (final_state, final_time, converged)
        - final_state: State vector at end of integration
        - final_time: Time at end of integration [days]
        - converged: True if steady state was reached
    """
    current_time = 0.0
    current_state = y0

    while current_time < max_time:
        # Integrate for check_interval
        t_span = (current_time, current_time + check_interval)
        result = self.solve(fun, t_span, current_state)

        # Get final state
        new_state = result.y[:, -1]

        # Check if steady state reached
        state_change = np.linalg.norm(new_state - current_state)
        state_norm = np.linalg.norm(new_state)

        if state_norm > 0:
            relative_change = state_change / state_norm
            if relative_change < steady_state_tol:
                return list(new_state), current_time + check_interval, True

        # Update for next iteration
        current_state = new_state
        current_time += check_interval

    # Max time exceeded without reaching steady state
    return list(current_state), current_time, False

pyadm1.core.solver.AdaptiveODESolver

Bases: ODESolver

Adaptive ODE solver that adjusts tolerances based on solution behavior.

Monitors the solution and can tighten tolerances if instabilities are detected, or relax them for faster computation when solution is smooth.

Source code in pyadm1/core/solver.py

class AdaptiveODESolver(ODESolver):
    """
    Adaptive ODE solver that adjusts tolerances based on solution behavior.

    Monitors the solution and can tighten tolerances if instabilities are
    detected, or relax them for faster computation when solution is smooth.
    """

    def __init__(
        self,
        config: Optional[SolverConfig] = None,
        adaptive: bool = True,
        min_rtol: float = 1e-8,
        max_rtol: float = 1e-4,
    ):
        """
        Initialize adaptive ODE solver.

        Args:
            config: Base solver configuration
            adaptive: Enable adaptive tolerance adjustment
            min_rtol: Minimum relative tolerance
            max_rtol: Maximum relative tolerance
        """
        super().__init__(config)
        self.adaptive = adaptive
        self.min_rtol = min_rtol
        self.max_rtol = max_rtol
        self._solution_history = []

    def solve(
        self,
        fun: Callable[[float, List[float]], List[float]],
        t_span: Tuple[float, float],
        y0: List[float],
        t_eval: Optional[np.ndarray] = None,
        dense_output: bool = False,
    ) -> "OdeResult":
        """
        Solve with adaptive tolerance adjustment.

        Same interface as parent class but monitors solution quality.
        """
        result = super().solve(fun, t_span, y0, t_eval, dense_output)

        if self.adaptive:
            self._update_tolerances(result)

        return result

    def _update_tolerances(self, result: "OdeResult") -> None:
        """
        Update solver tolerances based on solution behavior.

        Args:
            result: ODE solution result to analyze
        """
        # Calculate solution smoothness (using second derivatives)
        if len(result.t) > 2:
            y = result.y
            t = result.t

            # Estimate second derivatives
            dy_dt = np.gradient(y, t, axis=1)
            d2y_dt2 = np.gradient(dy_dt, t, axis=1)

            # Calculate maximum relative second derivative
            max_curvature = np.max(np.abs(d2y_dt2) / (np.abs(y) + 1e-10))

            # Adjust tolerances based on curvature
            if max_curvature > 1.0:
                # Solution has high curvature, tighten tolerances
                self.config.rtol = max(self.min_rtol, self.config.rtol * 0.5)
                self.config.atol = max(self.config.atol * 0.5, 1e-10)
            elif max_curvature < 0.1 and self.config.rtol < self.max_rtol:
                # Solution is smooth, can relax tolerances
                self.config.rtol = min(self.max_rtol, self.config.rtol * 1.5)
                self.config.atol = min(self.config.atol * 1.5, 1e-6)

Functions

__init__(config=None, adaptive=True, min_rtol=1e-08, max_rtol=0.0001)

Initialize adaptive ODE solver.

Parameters:

Name	Type	Description	Default
config	Optional[SolverConfig]	Base solver configuration	None
adaptive	bool	Enable adaptive tolerance adjustment	True
min_rtol	float	Minimum relative tolerance	1e-08
max_rtol	float	Maximum relative tolerance	0.0001

Source code in pyadm1/core/solver.py

def __init__(
    self,
    config: Optional[SolverConfig] = None,
    adaptive: bool = True,
    min_rtol: float = 1e-8,
    max_rtol: float = 1e-4,
):
    """
    Initialize adaptive ODE solver.

    Args:
        config: Base solver configuration
        adaptive: Enable adaptive tolerance adjustment
        min_rtol: Minimum relative tolerance
        max_rtol: Maximum relative tolerance
    """
    super().__init__(config)
    self.adaptive = adaptive
    self.min_rtol = min_rtol
    self.max_rtol = max_rtol
    self._solution_history = []

solve(fun, t_span, y0, t_eval=None, dense_output=False)

Solve with adaptive tolerance adjustment.

Same interface as parent class but monitors solution quality.

Source code in pyadm1/core/solver.py

def solve(
    self,
    fun: Callable[[float, List[float]], List[float]],
    t_span: Tuple[float, float],
    y0: List[float],
    t_eval: Optional[np.ndarray] = None,
    dense_output: bool = False,
) -> "OdeResult":
    """
    Solve with adaptive tolerance adjustment.

    Same interface as parent class but monitors solution quality.
    """
    result = super().solve(fun, t_span, y0, t_eval, dense_output)

    if self.adaptive:
        self._update_tolerances(result)

    return result

pyadm1.core.solver.SolverConfig dataclass

Configuration for ODE solver.

Attributes:

Name	Type	Description
method	str	Integration method ('BDF' for stiff ODEs)
rtol	float	Relative tolerance for solver
atol	float	Absolute tolerance for solver
min_step	float	Minimum allowed time step [days]
max_step	float	Maximum allowed time step [days]
first_step	Optional[float]	Initial step size [days]

Source code in pyadm1/core/solver.py

@dataclass
class SolverConfig:
    """
    Configuration for ODE solver.

    Attributes:
        method: Integration method ('BDF' for stiff ODEs)
        rtol: Relative tolerance for solver
        atol: Absolute tolerance for solver
        min_step: Minimum allowed time step [days]
        max_step: Maximum allowed time step [days]
        first_step: Initial step size [days]
    """

    method: str = "BDF"
    rtol: float = 1e-6
    atol: float = 1e-8
    min_step: float = 1e-6
    max_step: float = 0.1
    first_step: Optional[float] = None

pyadm1.core.solver.create_solver(method='BDF', adaptive=False, **kwargs)

Factory function to create appropriate solver instance.

Parameters:

Name	Type	Description	Default
method	str	Integration method ('BDF' recommended for ADM1)	'BDF'
adaptive	bool	Use adaptive tolerance adjustment	False
**kwargs		Additional configuration parameters for SolverConfig	{}

Returns:

Type	Description
ODESolver	Configured solver instance

Example

solver = create_solver(method='BDF', rtol=1e-7)

Use solver for integration...

Source code in pyadm1/core/solver.py

def create_solver(method: str = "BDF", adaptive: bool = False, **kwargs) -> ODESolver:
    """
    Factory function to create appropriate solver instance.

    Args:
        method: Integration method ('BDF' recommended for ADM1)
        adaptive: Use adaptive tolerance adjustment
        **kwargs: Additional configuration parameters for SolverConfig

    Returns:
        Configured solver instance

    Example:
        >>> solver = create_solver(method='BDF', rtol=1e-7)
        >>> # Use solver for integration...
    """
    config = SolverConfig(method=method, **kwargs)

    if adaptive:
        return AdaptiveODESolver(config)
    else:
        return ODESolver(config)