Core ADM1 Model

pyadm1.core.adm1.ADM1

Main class implementing ADM1 as pure ODE system.

This class manages the ADM1 state, parameters, and provides methods for simulation including influent stream creation, gas production calculation, and state tracking.

Attributes:

Name	Type	Description
V_liq		Liquid volume [m³]
T_ad	float	Operating temperature [K]
feedstock	Feedstock	Feedstock object for substrate management

Example

feedstock = Feedstock(feeding_freq=48) adm1 = ADM1(feedstock, V_liq=2000, T_ad=308.15) adm1.create_influent([15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 0)

Source code in pyadm1/core/adm1.py

class ADM1:
    """
    Main class implementing ADM1 as pure ODE system.

    This class manages the ADM1 state, parameters, and provides methods for
    simulation including influent stream creation, gas production calculation,
    and state tracking.

    Attributes:
        V_liq: Liquid volume [m³]
        T_ad: Operating temperature [K]
        feedstock: Feedstock object for substrate management

    Example:
        >>> feedstock = Feedstock(feeding_freq=48)
        >>> adm1 = ADM1(feedstock, V_liq=2000, T_ad=308.15)
        >>> adm1.create_influent([15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 0)
    """

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

        Args:
            feedstock: Feedstock object containing substrate information. E.g. used to calculate ADM1 input stream.
            V_liq: Liquid volume [m³]
            V_gas: Gas volume [m³]
            T_ad: Operating temperature [K] (default: 308.15 = 35°C)
        """
        # Physical parameters
        self.V_liq = V_liq  # liquid volume of digester
        self._V_gas = V_gas  # gas volume of digester
        self._V_ad = self.V_liq + self._V_gas  # total volume of digester: liquid + gas volume
        self._T_ad = T_ad  # temperature inside the digester

        # Constants
        self._R = 0.08314  # 0.083145  Gas constant [bar·M^-1·K^-1]
        self._T_base = 298.15  # Base temperature [K] 25°C
        self._p_atm = 1.013  # Atmospheric pressure [bar]

        # Calculated parameters
        self._RT = self._R * self._T_ad  # R * T_ad
        # external pressures
        self._p_ext = self._p_atm - 0.0084147 * np.exp(0.054 * (self._T_ad - 273.15))

        # Feedstock and state
        self._feedstock = feedstock
        # vector of volumetric flow rates of the substrates. Length must be equal to the number of
        # substrates defined in xml
        self._Q: Optional[List[float]] = None
        self._state_input: Optional[List[float]] = None  # contains ADM1 input stream as a 34dim vector

        # Result tracking lists
        # Result tracking lists
        self._Q_GAS: List[float] = []  # produced biogas over all simulations [m³/d]
        self._Q_CH4: List[float] = []  # produced methane over all simulations [m³/d]
        self._Q_CO2: List[float] = []  # produced CO2 over all simulations [m³/d]
        self._P_GAS: List[float] = []  # gas pressures over all simulations [bar]
        self._pH_l: List[float] = []  # pH values over all simulations [-]
        self._FOSTAC: List[float] = []  # ratio of VFA over TA over all simulations [-]
        self._AcvsPro: List[float] = []  # ratio of acetic over propionic acid over all simulations [-]
        self._VFA: List[float] = []  # VFA concentrations over all simulations [g/L]
        self._TAC: List[float] = []  # TA concentrations over all simulations [g CaCO3/L]

    def create_influent(self, Q: List[float], i: int) -> None:
        """
        Create ADM1 input stream from volumetric flow rates.

        Calculates the ADM1 influent state by mixing substrate streams according
        to their volumetric flow rates. The resulting influent composition is
        stored internally for use in ODE calculations.

        Args:
            Q: Volumetric flow rates for each substrate [m³/d]
                Length must equal number of substrates in feedstock
            i: Time step index for accessing influent dataframe

        Example:
            >>> adm1.create_influent([15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 0)
        """
        self._Q = Q
        influent_state = self._feedstock.get_influent_dataframe(Q)
        self._set_influent(influent_state, i)

    def calc_gas(self, pi_Sh2: float, pi_Sch4: float, pi_Sco2: float, pTOTAL: float) -> Tuple[float, float, float, float]:
        """
        Calculate biogas production rates from partial pressures.

        Uses the ideal gas law and Henry's constants to calculate gas flow rates
        from the gas phase partial pressures.

        Args:
            pi_Sh2: Hydrogen partial pressure [bar]
            pi_Sch4: Methane partial pressure [bar]
            pi_Sco2: CO2 partial pressure [bar]
            pTOTAL: Total gas pressure [bar]

        Returns:
            Tuple containing:
                - q_gas: Total biogas flow rate [m³/d]
                - q_ch4: Methane flow rate [m³/d]
                - q_co2: CO2 flow rate [m³/d]
                - p_gas: Total gas partial pressure (excl. H2O) [bar]

        Example:
            >>> q_gas, q_ch4, q_co2, p_gas = adm1.calc_gas(5e-6, 0.55, 0.42, 0.98)
            >>> print(f"Biogas: {q_gas:.1f} m³/d, Methane: {q_ch4:.1f} m³/d")
        """
        _, k_p, _, _, _, _ = ADMParams.getADMgasparams(self._R, self._T_base, self._T_ad)

        # Ideal gas law constant for conversion
        NQ = 44.643

        # Total biogas flow from pressure difference
        q_gas = k_p * (pTOTAL - self._p_ext) / (self._RT / 1000 * NQ) * self.V_liq

        # Total gas partial pressure (excluding water vapor)
        p_gas = pi_Sh2 + pi_Sch4 + pi_Sco2

        # Ensure non-negative gas flows
        if isinstance(q_gas, np.ndarray):
            q_gas = np.maximum(q_gas, 0.0)
        else:
            q_gas = max(q_gas, 0.0)

        # Calculate component flows from partial pressures
        if p_gas > 0:
            q_ch4 = q_gas * (pi_Sch4 / p_gas)
            q_co2 = q_gas * (pi_Sco2 / p_gas)
        else:
            q_ch4 = 0.0
            q_co2 = 0.0

        # Ensure non-negative
        if isinstance(q_ch4, np.ndarray):
            q_ch4 = np.maximum(q_ch4, 0.0)
            q_co2 = np.maximum(q_co2, 0.0)
        else:
            q_ch4 = max(q_ch4, 0.0)
            q_co2 = max(q_co2, 0.0)

        return q_gas, q_ch4, q_co2, p_gas

    def resume_from_broken_simulation(self, Q_CH4):
        for Qch4 in Q_CH4:
            self._Q_CH4.append(Qch4)

    def save_final_state_in_csv(self, simulate_results: List[List[float]], filename: str = "digester_final.csv") -> None:
        """
        Save final ADM1 state vector to CSV file.

        Exports only the last state from simulation results, which can be used
        as initial state for subsequent simulations.

        Args:
            simulate_results: List of ADM1 state vectors from simulation
            filename: Output CSV filename

        Example:
            >>> results = [[0.01]*37, [0.02]*37, [0.03]*37]
            >>> adm1.save_final_state_in_csv(results, 'final_state.csv')
        """
        columns = [
            *self._feedstock.header()[:-1],
            "pi_Sh2",
            "pi_Sch4",
            "pi_Sco2",
            "pTOTAL",
        ]

        simulate_results_df = pd.DataFrame.from_records(simulate_results)
        simulate_results_df.columns = columns

        # Keep only the last row
        last_simulated_result = simulate_results_df[-1:]
        last_simulated_result.to_csv(filename, index=False)

    def print_params_at_current_state(self, state_ADM1xp: List[float]) -> None:
        """
        Calculate and store process parameters from current state.

        Computes key process indicators including pH, VFA, TAC,
        and gas production rates. Also stores values in tracking lists.

        Args:
            state_ADM1xp: Current ADM1 state vector (37 elements)

        Example:
            >>> adm1.print_params_at_current_state(state_vector)
        """
        # Calculate process indicators using DLL
        self._pH_l.append(np.round(ADMstate.calcPHOfADMstate(state_ADM1xp), 1))
        self._FOSTAC.append(np.round(ADMstate.calcFOSTACOfADMstate(state_ADM1xp).Value, 2))
        self._AcvsPro.append(np.round(ADMstate.calcAcetic_vs_PropionicOfADMstate(state_ADM1xp).Value, 1))
        self._VFA.append(np.round(ADMstate.calcVFAOfADMstate(state_ADM1xp, "gHAceq/l").Value, 2))
        self._TAC.append(np.round(ADMstate.calcTACOfADMstate(state_ADM1xp, "gCaCO3eq/l").Value, 1))

        # Ensure at least 2 values in lists (because the last three values go to the controller)
        # I am assuming here that we start from a steady state
        if len(self._pH_l) < 2:
            self._pH_l.append(self._pH_l[-1])
            self._FOSTAC.append(self._FOSTAC[-1])
            self._AcvsPro.append(self._AcvsPro[-1])
            self._VFA.append(self._VFA[-1])
            self._TAC.append(self._TAC[-1])

        # Calculate and store gas production
        q_gas, q_ch4, q_co2, p_gas = self.calc_gas(state_ADM1xp[33], state_ADM1xp[34], state_ADM1xp[35], state_ADM1xp[36])

        self._Q_GAS.append(q_gas)
        self._Q_CH4.append(q_ch4)
        self._Q_CO2.append(q_co2)
        self._P_GAS.append(p_gas)

        # Ensure at least 2 values, because the last three values go to controller.
        # I am assuming here that we start from a steady state
        if len(self._Q_GAS) < 2:
            for _ in range(2):  # to have at least 4 values in the lists
                self._Q_GAS.append(q_gas)
                self._Q_CH4.append(q_ch4)
                self._Q_CO2.append(q_co2)
                self._P_GAS.append(p_gas)

    def ADM1_ODE(self, t: float, state_zero: List[float]) -> Tuple[float, ...]:
        """
        Calculate derivatives for ADM1 ODE system.

        This is the main ODE function that computes dy/dt for all 37 state
        variables. Uses process rate equations and stoichiometric relationships.

        Args:
            t: Current time [days] (not used, system is autonomous)
            state_zero: Current ADM1 state vector (37 elements)

        Returns:
            Tuple of 37 derivatives (dy/dt)

        Note:
            This method is called by the ODE solver and should not be called
            directly by users.
        """
        # Get all ADM1 parameters
        params_tuple = ADMParams.getADMparams(self._R, self._T_base, self._T_ad)

        # Get substrate-dependent parameters
        substrate_params = self._get_substrate_dependent_params()

        # Unpack frequently used parameters
        (
            N_xc,
            N_I,
            N_aa,
            C_xc,
            C_sI,
            C_ch,
            C_pr,
            C_li,
            C_xI,
            C_su,
            C_aa,
            f_fa_li,
            C_fa,
            f_h2_su,
            f_bu_su,
            f_pro_su,
            f_ac_su,
            N_bac,
            C_bu,
            C_pro,
            C_ac,
            C_bac,
            Y_su,
            f_h2_aa,
            f_va_aa,
            f_bu_aa,
            f_pro_aa,
            f_ac_aa,
            C_va,
            Y_aa,
            Y_fa,
            Y_c4,
            Y_pro,
            C_ch4,
            Y_ac,
            Y_h2,
        ) = params_tuple[:36]

        # Get kinetic parameters (starting from index 36)
        kinetic_params = {
            "K_S_IN": params_tuple[36],
            "k_m_su": params_tuple[37],
            "K_S_su": params_tuple[38],
            "k_m_aa": params_tuple[41],
            "K_S_aa": params_tuple[42],
            "k_m_fa": params_tuple[43],
            "K_S_fa": params_tuple[44],
            "K_I_h2_fa": params_tuple[45],
            "k_m_c4": params_tuple[46],
            "K_S_c4": params_tuple[47],
            "K_I_h2_c4": params_tuple[48],
            "k_m_pro": params_tuple[49],
            "K_S_pro": params_tuple[50],
            "K_I_h2_pro": params_tuple[51],
            "k_m_ac": params_tuple[52],
            "K_S_ac": params_tuple[53],
            "K_I_nh3": params_tuple[54],
            "k_m_h2": params_tuple[57],
            "K_S_h2": params_tuple[58],
            "k_dec_X_su": params_tuple[61],
            "k_dec_X_aa": params_tuple[62],
            "k_dec_X_fa": params_tuple[63],
            "k_dec_X_c4": params_tuple[64],
            "k_dec_X_pro": params_tuple[65],
            "k_dec_X_ac": params_tuple[66],
            "k_dec_X_h2": params_tuple[67],
        }

        # Get acid-base and gas parameters
        acid_base_params = {
            "K_w": params_tuple[68],
            "K_a_va": params_tuple[69],
            "K_a_bu": params_tuple[70],
            "K_a_pro": params_tuple[71],
            "K_a_ac": params_tuple[72],
            "K_a_co2": params_tuple[73],
            "K_a_IN": params_tuple[74],
            "k_A_B_va": params_tuple[75],
            "k_A_B_bu": params_tuple[76],
            "k_A_B_pro": params_tuple[77],
            "k_A_B_ac": params_tuple[78],
            "k_A_B_co2": params_tuple[79],
            "k_A_B_IN": params_tuple[80],
        }

        gas_params = {
            "k_p": params_tuple[82],
            "k_L_a": params_tuple[83],
            "K_H_co2": params_tuple[84],
            "K_H_ch4": params_tuple[85],
            "K_H_h2": params_tuple[86],
        }

        # Get pH and inhibition parameters
        K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2 = ADMParams.getADMinhibitionparams()

        # Calculate pH and H+ concentration
        pH = ADMstate.calcPHOfADMstate(state_zero)
        S_H_ion = 10 ** (-pH)

        # Unpack state variables
        S_su, S_aa, S_fa, S_va, S_bu, S_pro, S_ac, S_h2 = state_zero[0:8]
        S_ch4, S_co2, S_nh4_ion, S_I = state_zero[8:12]
        X_xc, X_ch, X_pr, X_li = state_zero[12:16]
        X_su, X_aa, X_fa, X_c4, X_pro, X_ac, X_h2, X_I, X_p = state_zero[16:25]
        S_cation, S_anion = state_zero[25:27]
        S_va_ion, S_bu_ion, S_pro_ion, S_ac_ion = state_zero[27:31]
        S_hco3_ion, S_nh3 = state_zero[31:33]
        p_gas_h2, p_gas_ch4, p_gas_co2, pTOTAL = state_zero[33:37]

        # Get influent concentrations
        q_ad = np.sum(self._Q) if self._Q is not None else 0.0
        state_in = self._state_input if self._state_input is not None else [0.0] * 34

        # Calculate inhibition factors
        bio = BiochemicalProcesses()
        inhibitions = bio.calculate_inhibition_factors(
            S_H_ion,
            S_h2,
            S_nh4_ion,
            S_nh3,
            K_pH_aa,
            nn_aa,
            K_pH_ac,
            n_ac,
            K_pH_h2,
            n_h2,
            kinetic_params["K_S_IN"],
            kinetic_params["K_I_h2_fa"],
            kinetic_params["K_I_h2_c4"],
            kinetic_params["K_I_h2_pro"],
            kinetic_params["K_I_nh3"],
        )

        # Calculate biochemical process rates
        process_rates = bio.calculate_process_rates(state_zero, inhibitions, kinetic_params, substrate_params)
        Rho_1, Rho_2, Rho_3, Rho_4, Rho_5, Rho_6, Rho_7, Rho_8, Rho_9 = process_rates[:9]
        Rho_10, Rho_11, Rho_12, Rho_13, Rho_14, Rho_15, Rho_16 = process_rates[9:16]
        Rho_17, Rho_18, Rho_19 = process_rates[16:19]

        # Calculate acid-base rates
        acid_base_rates = bio.calculate_acid_base_rates(state_zero, acid_base_params)
        Rho_A_4, Rho_A_5, Rho_A_6, Rho_A_7, Rho_A_10, Rho_A_11 = acid_base_rates

        # Calculate gas transfer rates
        gas_rates = bio.calculate_gas_transfer_rates(state_zero, gas_params, self._RT, self.V_liq, self._V_gas, self._p_ext)
        Rho_T_8, Rho_T_9, Rho_T_10, Rho_T_11 = gas_rates

        # Extract substrate-dependent fractions
        f_ch_xc, f_pr_xc, f_li_xc, f_xI_xc, f_sI_xc, f_xp_xc, _, _, _, _, _, _, _, _ = substrate_params.values()

        # Calculate biomass decay fractions
        f_p = 0.08
        f_ch_xb = f_ch_xc / (f_ch_xc + f_pr_xc + f_li_xc) * (1 - f_p)
        f_pr_xb = f_pr_xc / (f_ch_xc + f_pr_xc + f_li_xc) * (1 - f_p)
        f_li_xb = f_li_xc / (f_ch_xc + f_pr_xc + f_li_xc) * (1 - f_p)

        # Differential equations for soluble components (1-12)
        diff_S_su = q_ad / self.V_liq * (state_in[0] - S_su) + Rho_2 + (1 - f_fa_li) * Rho_4 - Rho_5

        diff_S_aa = q_ad / self.V_liq * (state_in[1] - S_aa) + Rho_3 - Rho_6

        diff_S_fa = q_ad / self.V_liq * (state_in[2] - S_fa) + f_fa_li * Rho_4 - Rho_7

        diff_S_va = q_ad / self.V_liq * (state_in[3] - S_va) + (1 - Y_aa) * f_va_aa * Rho_6 - Rho_8

        diff_S_bu = (
            q_ad / self.V_liq * (state_in[4] - S_bu) + (1 - Y_su) * f_bu_su * Rho_5 + (1 - Y_aa) * f_bu_aa * Rho_6 - Rho_9
        )

        diff_S_pro = (
            q_ad / self.V_liq * (state_in[5] - S_pro)
            + (1 - Y_su) * f_pro_su * Rho_5
            + (1 - Y_aa) * f_pro_aa * Rho_6
            + (1 - Y_c4) * 0.54 * Rho_8
            - Rho_10
        )

        diff_S_ac = (
            q_ad / self.V_liq * (state_in[6] - S_ac)
            + (1 - Y_su) * f_ac_su * Rho_5
            + (1 - Y_aa) * f_ac_aa * Rho_6
            + (1 - Y_fa) * 0.7 * Rho_7
            + (1 - Y_c4) * 0.31 * Rho_8
            + (1 - Y_c4) * 0.8 * Rho_9
            + (1 - Y_pro) * 0.57 * Rho_10
            - Rho_11
        )

        diff_S_h2 = (
            q_ad / self.V_liq * (state_in[7] - S_h2)
            + (1 - Y_su) * f_h2_su * Rho_5
            + (1 - Y_aa) * f_h2_aa * Rho_6
            + (1 - Y_fa) * 0.3 * Rho_7
            + (1 - Y_c4) * 0.15 * Rho_8
            + (1 - Y_c4) * 0.2 * Rho_9
            + (1 - Y_pro) * 0.43 * Rho_10
            - Rho_12
            - self._V_gas / self.V_liq * Rho_T_8
        )

        diff_S_ch4 = (
            q_ad / self.V_liq * (state_in[8] - S_ch4)
            + (1 - Y_ac) * Rho_11
            + (1 - Y_h2) * Rho_12
            - self._V_gas / self.V_liq * Rho_T_9
        )

        # CO2 balance with carbon stoichiometry
        s_1 = -C_xc + f_sI_xc * C_sI + f_ch_xc * C_ch + f_pr_xc * C_pr + f_li_xc * C_li + f_xI_xc * C_xI
        s_2 = -C_ch + C_su
        s_3 = -C_pr + C_aa
        s_4 = -C_li + (1 - f_fa_li) * C_su + f_fa_li * C_fa
        s_5 = -C_su + (1 - Y_su) * (f_bu_su * C_bu + f_pro_su * C_pro + f_ac_su * C_ac) + Y_su * C_bac
        s_6 = -C_aa + (1 - Y_aa) * (f_va_aa * C_va + f_bu_aa * C_bu + f_pro_aa * C_pro + f_ac_aa * C_ac) + Y_aa * C_bac
        s_7 = -C_fa + (1 - Y_fa) * 0.7 * C_ac + Y_fa * C_bac
        s_8 = -C_va + (1 - Y_c4) * 0.54 * C_pro + (1 - Y_c4) * 0.31 * C_ac + Y_c4 * C_bac
        s_9 = -C_bu + (1 - Y_c4) * 0.8 * C_ac + Y_c4 * C_bac
        s_10 = -C_pro + (1 - Y_pro) * 0.57 * C_ac + Y_pro * C_bac
        s_11 = -C_ac + (1 - Y_ac) * C_ch4 + Y_ac * C_bac
        s_12 = (1 - Y_h2) * C_ch4 + Y_h2 * C_bac
        s_13 = -C_bac + C_xc

        Sigma = (
            s_1 * Rho_1
            + s_2 * Rho_2
            + s_3 * Rho_3
            + s_4 * Rho_4
            + s_5 * Rho_5
            + s_6 * Rho_6
            + s_7 * Rho_7
            + s_8 * Rho_8
            + s_9 * Rho_9
            + s_10 * Rho_10
            + s_11 * Rho_11
            + s_12 * Rho_12
            + s_13 * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
        )

        diff_S_co2 = q_ad / self.V_liq * (state_in[9] - S_co2) - Sigma - self._V_gas / self.V_liq * Rho_T_10 + Rho_A_10

        # Nitrogen balance
        diff_S_nh4_ion = (
            q_ad / self.V_liq * (state_in[10] - S_nh4_ion)
            - Y_su * N_bac * Rho_5
            + (N_aa - Y_aa * N_bac) * Rho_6
            - Y_fa * N_bac * Rho_7
            - Y_c4 * N_bac * Rho_8
            - Y_c4 * N_bac * Rho_9
            - Y_pro * N_bac * Rho_10
            - Y_ac * N_bac * Rho_11
            - Y_h2 * N_bac * Rho_12
            + (N_bac - N_xc) * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
            + Rho_A_11
        )

        diff_S_I = q_ad / self.V_liq * (state_in[11] - S_I) + f_sI_xc * Rho_1

        # Differential equations for particulate components (13-24)
        diff_X_xc = q_ad / self.V_liq * (state_in[12] - X_xc) - Rho_1

        diff_X_ch = (
            q_ad / self.V_liq * (state_in[13] - X_ch)
            + f_ch_xc * Rho_1
            - Rho_2
            + f_ch_xb * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
        )

        diff_X_pr = (
            q_ad / self.V_liq * (state_in[14] - X_pr)
            + f_pr_xc * Rho_1
            - Rho_3
            + f_pr_xb * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
        )

        diff_X_li = (
            q_ad / self.V_liq * (state_in[15] - X_li)
            + f_li_xc * Rho_1
            - Rho_4
            + f_li_xb * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
        )

        diff_X_su = q_ad / self.V_liq * (state_in[16] - X_su) + Y_su * Rho_5 - Rho_13
        diff_X_aa = q_ad / self.V_liq * (state_in[17] - X_aa) + Y_aa * Rho_6 - Rho_14
        diff_X_fa = q_ad / self.V_liq * (state_in[18] - X_fa) + Y_fa * Rho_7 - Rho_15
        diff_X_c4 = q_ad / self.V_liq * (state_in[19] - X_c4) + Y_c4 * Rho_8 + Y_c4 * Rho_9 - Rho_16
        diff_X_pro = q_ad / self.V_liq * (state_in[20] - X_pro) + Y_pro * Rho_10 - Rho_17
        diff_X_ac = q_ad / self.V_liq * (state_in[21] - X_ac) + Y_ac * Rho_11 - Rho_18
        diff_X_h2 = q_ad / self.V_liq * (state_in[22] - X_h2) + Y_h2 * Rho_12 - Rho_19
        diff_X_I = q_ad / self.V_liq * (state_in[23] - X_I) + f_xI_xc * Rho_1

        diff_X_p = (
            q_ad / self.V_liq * (state_in[24] - X_p)
            - f_xp_xc * Rho_1
            + f_p * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
        )

        # Differential equations for ions (25-32)
        diff_S_cation = q_ad / self.V_liq * (state_in[25] - S_cation)
        diff_S_anion = q_ad / self.V_liq * (state_in[26] - S_anion)
        diff_S_va_ion = q_ad / self.V_liq * (state_in[27] - S_va_ion) - Rho_A_4
        diff_S_bu_ion = q_ad / self.V_liq * (state_in[28] - S_bu_ion) - Rho_A_5
        diff_S_pro_ion = q_ad / self.V_liq * (state_in[29] - S_pro_ion) - Rho_A_6
        diff_S_ac_ion = q_ad / self.V_liq * (state_in[30] - S_ac_ion) - Rho_A_7
        diff_S_hco3_ion = q_ad / self.V_liq * (state_in[31] - S_hco3_ion) - Rho_A_10
        diff_S_nh3 = q_ad / self.V_liq * (state_in[32] - S_nh3) - Rho_A_11

        # Differential equations for gas phase (33-36)
        diff_p_gas_h2 = Rho_T_8 * self._RT / 16 - p_gas_h2 / pTOTAL * Rho_T_11
        diff_p_gas_ch4 = Rho_T_9 * self._RT / 64 - p_gas_ch4 / pTOTAL * Rho_T_11
        diff_p_gas_co2 = Rho_T_10 * self._RT - p_gas_co2 / pTOTAL * Rho_T_11
        diff_pTOTAL = self._RT / 16 * Rho_T_8 + self._RT / 64 * Rho_T_9 + self._RT * Rho_T_10 - 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_ion,
            diff_S_I,
            diff_X_xc,
            diff_X_ch,
            diff_X_pr,
            diff_X_li,
            diff_X_su,
            diff_X_aa,
            diff_X_fa,
            diff_X_c4,
            diff_X_pro,
            diff_X_ac,
            diff_X_h2,
            diff_X_I,
            diff_X_p,
            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_ion,
            diff_S_nh3,
            diff_p_gas_h2,
            diff_p_gas_ch4,
            diff_p_gas_co2,
            diff_pTOTAL,
        )

    def _set_influent(self, influent_state: pd.DataFrame, i: int) -> None:
        """
        Set influent values from dataframe.

        Internal method to extract influent state at time index i and store
        it for use in ODE calculations.

        Args:
            influent_state: DataFrame with ADM1 input variables over time
            i: Time step index (uses last row if i exceeds available data)
        """
        # Handle index out of bounds by using last row (steady-state assumption)
        max_index = len(influent_state) - 1
        if i > max_index:
            i = max_index

        # Extract all state variables at time index i
        self._state_input = [
            influent_state["S_su"].iloc[i],  # kg COD.m^-3
            influent_state["S_aa"].iloc[i],  # kg COD.m^-3
            influent_state["S_fa"].iloc[i],  # kg COD.m^-3
            influent_state["S_va"].iloc[i],  # kg COD.m^-3
            influent_state["S_bu"].iloc[i],  # kg COD.m^-3
            influent_state["S_pro"].iloc[i],  # kg COD.m^-3
            influent_state["S_ac"].iloc[i],  # kg COD.m^-3
            influent_state["S_h2"].iloc[i],  # kg COD.m^-3
            influent_state["S_ch4"].iloc[i],  # kg COD.m^-3
            influent_state["S_co2"].iloc[i],  # kmole C.m^-3 (S_IC_in)
            influent_state["S_nh4"].iloc[i],  # kmole N.m^-3 (S_IN_in)
            influent_state["S_I"].iloc[i],  # kg COD.m^-3
            influent_state["X_xc"].iloc[i],  # kg COD.m^-3
            influent_state["X_ch"].iloc[i],  # kg COD.m^-3
            influent_state["X_pr"].iloc[i],  # kg COD.m^-3
            influent_state["X_li"].iloc[i],  # kg COD.m^-3
            influent_state["X_su"].iloc[i],  # kg COD.m^-3
            influent_state["X_aa"].iloc[i],  # kg COD.m^-3
            influent_state["X_fa"].iloc[i],  # kg COD.m^-3
            influent_state["X_c4"].iloc[i],  # kg COD.m^-3
            influent_state["X_pro"].iloc[i],  # kg COD.m^-3
            influent_state["X_ac"].iloc[i],  # kg COD.m^-3
            influent_state["X_h2"].iloc[i],  # kg COD.m^-3
            influent_state["X_I"].iloc[i],  # kg COD.m^-3
            influent_state["X_p"].iloc[i],  # kg COD.m^-3
            influent_state["S_cation"].iloc[i],  # kmole.m^-3
            influent_state["S_anion"].iloc[i],  # kmole.m^-3
            influent_state["S_va_ion"].iloc[i],  # kg COD.m^-3
            influent_state["S_bu_ion"].iloc[i],  # kg COD.m^-3
            influent_state["S_pro_ion"].iloc[i],  # kg COD.m^-3
            influent_state["S_ac_ion"].iloc[i],  # kg COD.m^-3
            influent_state["S_hco3_ion"].iloc[i],  # kg COD.m^-3
            influent_state["S_nh3"].iloc[i],  # kg COD.m^-3
            influent_state["Q"].iloc[i],  # m^3/d (q_ad)
        ]

    def _get_substrate_dependent_params(self) -> dict:
        """
        Get substrate-dependent parameters from C# DLL with calibration override support.

        Calculates ADM1 parameters that depend on substrate composition using weighted
        averaging based on volumetric flow rates. If calibration parameters are set
        (via _calibration_params attribute), those values override the calculated ones.
        TODO: do I want this behaviour that during calibration the parameters are overriden?

        Returns:
            dict: Substrate-dependent parameters with keys:
                - f_ch_xc, f_pr_xc, f_li_xc: Composite fractions
                - f_xI_xc, f_sI_xc, f_xp_xc: Inert and product fractions
                - k_dis: Disintegration rate [1/d]
                - k_hyd_ch, k_hyd_pr, k_hyd_li: Hydrolysis rates [1/d]
                - k_m_c4, k_m_pro, k_m_ac, k_m_h2: Uptake rates [1/d]

        Example:
            >>> # Normal usage
            >>> params = adm1._get_substrate_dependent_params()
            >>>
            >>> # With calibration override
            >>> adm1._calibration_params = {'k_dis': 0.55, 'Y_su': 0.105}
            >>> params = adm1._get_substrate_dependent_params()
            >>> # k_dis will be 0.55 instead of calculated value
        """
        if self._Q is None:
            # Return default values if Q not set
            base_params = {
                "f_ch_xc": 0.2,
                "f_pr_xc": 0.2,
                "f_li_xc": 0.3,
                "f_xI_xc": 0.2,
                "f_sI_xc": 0.1,
                "f_xp_xc": 0.0,
                "k_dis": 0.5,
                "k_hyd_ch": 10.0,
                "k_hyd_pr": 10.0,
                "k_hyd_li": 10.0,
                "k_m_c4": 20.0,
                "k_m_pro": 13.0,
                "k_m_ac": 8.0,
                "k_m_h2": 35.0,
            }
        else:

            def _calc_params_for_q(q_arg):
                f_ch_xc, f_pr_xc, f_li_xc, f_xI_xc, f_sI_xc, f_xp_xc = self._feedstock.mySubstrates().calcfFactors(q_arg)
                f_xp_xc = max(f_xp_xc, 0.0)

                k_dis = self._feedstock.mySubstrates().calcDisintegrationParam(q_arg)
                k_hyd_ch, k_hyd_pr, k_hyd_li = self._feedstock.mySubstrates().calcHydrolysisParams(q_arg)
                k_m_c4, k_m_pro, k_m_ac, k_m_h2 = self._feedstock.mySubstrates().calcMaxUptakeRateParams(q_arg)

                return {
                    "f_ch_xc": f_ch_xc,
                    "f_pr_xc": f_pr_xc,
                    "f_li_xc": f_li_xc,
                    "f_xI_xc": f_xI_xc,
                    "f_sI_xc": f_sI_xc,
                    "f_xp_xc": f_xp_xc,
                    "k_dis": k_dis,
                    "k_hyd_ch": k_hyd_ch,
                    "k_hyd_pr": k_hyd_pr,
                    "k_hyd_li": k_hyd_li,
                    "k_m_c4": k_m_c4,
                    "k_m_pro": k_m_pro,
                    "k_m_ac": k_m_ac,
                    "k_m_h2": k_m_h2,
                }

            errors = []

            # Preferred path for pythonnet 3.x on Linux/Mono: pass explicit System.Double[,]
            try:
                from System import Array, Double

                q_values = [float(q) for q in self._Q]
                q_2d = Array.CreateInstance(Double, 1, len(q_values))
                for idx, value in enumerate(q_values):
                    q_2d[0, idx] = value

                base_params = _calc_params_for_q(q_2d)
            except Exception as exc:
                errors.append(f"System.Double[,] path failed: {exc}")

                # Fallback path: numpy 2D array works on some local runtimes
                try:
                    q_2d_np = np.atleast_2d(np.asarray(self._Q, dtype=float))
                    base_params = _calc_params_for_q(q_2d_np)
                except Exception as exc_np:
                    errors.append(f"numpy 2D path failed: {exc_np}")

                    # Final fallback for legacy runtimes that still accept 1D arrays
                    try:
                        q_1d = [float(q) for q in self._Q]
                        base_params = _calc_params_for_q(q_1d)
                    except Exception as exc_1d:
                        errors.append(f"1D path failed: {exc_1d}")
                        raise TypeError("Failed to evaluate substrate-dependent parameters. " + " | ".join(errors))

        # Apply calibration parameter overrides if they exist
        if hasattr(self, "_calibration_params") and self._calibration_params:
            for param_name, param_value in self._calibration_params.items():
                if param_name in base_params:
                    base_params[param_name] = param_value

        return base_params

    def set_calibration_parameters(self, parameters: dict) -> None:
        """
        Set calibration parameters that override substrate-dependent calculations.

        Args:
            parameters: Parameter values as {param_name: value}.

        Example:
            >>> adm1.set_calibration_parameters({
            ...     'k_dis': 0.55,
            ...     'k_hyd_ch': 11.0,
            ...     'Y_su': 0.105
            ... })
        """
        if not hasattr(self, "_calibration_params"):
            self._calibration_params = {}

        self._calibration_params.update(parameters)

    def clear_calibration_parameters(self) -> None:
        """
        Clear all calibration parameters and revert to substrate-dependent calculations.

        Example:
            >>> adm1.clear_calibration_parameters()
        """
        if hasattr(self, "_calibration_params"):
            self._calibration_params = {}

    def get_calibration_parameters(self) -> dict:
        """
        Get currently set calibration parameters.

        Returns:
            dict: Current calibration parameters as {param_name: value}.

        Example:
            >>> params = adm1.get_calibration_parameters()
            >>> print(params)
            {'k_dis': 0.55, 'Y_su': 0.105}
        """
        if hasattr(self, "_calibration_params"):
            return self._calibration_params.copy()
        return {}

    # Properties for accessing model parameters and results
    @property
    def T_ad(self) -> float:
        """Operating temperature [K]."""
        return self._T_ad

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

    @property
    def Q_GAS(self) -> List[float]:
        """Biogas production rates over all simulations [m³/d]."""
        return self._Q_GAS

    @property
    def Q_CH4(self) -> List[float]:
        """Methane production rates over all simulations [m³/d]."""
        return self._Q_CH4

    @property
    def Q_CO2(self) -> List[float]:
        """CO2 production rates over all simulations [m³/d]."""
        return self._Q_CO2

    @property
    def P_GAS(self) -> List[float]:
        """Gas pressures over all simulations [bar]."""
        return self._P_GAS

    @property
    def pH_l(self) -> List[float]:
        """pH values over all simulations."""
        return self._pH_l

    @property
    def VFA_TA(self) -> List[float]:
        """VFA/TA ratios over all simulations."""
        return self._FOSTAC

    @property
    def AcvsPro(self) -> List[float]:
        """Acetic/Propionic acid ratios over all simulations."""
        return self._AcvsPro

    @property
    def VFA(self) -> List[float]:
        """VFA concentrations over all simulations [g/L]."""
        return self._VFA

    @property
    def TAC(self) -> List[float]:
        """TA concentrations over all simulations [g CaCO3 eq/L]."""
        return self._TAC

Attributes

AcvsPro property

Acetic/Propionic acid ratios over all simulations.

P_GAS property

Gas pressures over all simulations [bar].

Q_CH4 property

Methane production rates over all simulations [m³/d].

Q_CO2 property

CO2 production rates over all simulations [m³/d].

Q_GAS property

Biogas production rates over all simulations [m³/d].

TAC property

TA concentrations over all simulations [g CaCO3 eq/L].

T_ad property

Operating temperature [K].

VFA property

VFA concentrations over all simulations [g/L].

VFA_TA property

VFA/TA ratios over all simulations.

feedstock property

Feedstock object.

pH_l property

pH values over all simulations.

Functions

ADM1_ODE(t, state_zero)

Calculate derivatives for ADM1 ODE system.

This is the main ODE function that computes dy/dt for all 37 state variables. Uses process rate equations and stoichiometric relationships.

Parameters:

Name	Type	Description	Default
t	float	Current time [days] (not used, system is autonomous)	required
state_zero	List[float]	Current ADM1 state vector (37 elements)	required

Returns:

Type	Description
Tuple[float, ...]	Tuple of 37 derivatives (dy/dt)

Note

This method is called by the ODE solver and should not be called directly by users.

Source code in pyadm1/core/adm1.py

def ADM1_ODE(self, t: float, state_zero: List[float]) -> Tuple[float, ...]:
    """
    Calculate derivatives for ADM1 ODE system.

    This is the main ODE function that computes dy/dt for all 37 state
    variables. Uses process rate equations and stoichiometric relationships.

    Args:
        t: Current time [days] (not used, system is autonomous)
        state_zero: Current ADM1 state vector (37 elements)

    Returns:
        Tuple of 37 derivatives (dy/dt)

    Note:
        This method is called by the ODE solver and should not be called
        directly by users.
    """
    # Get all ADM1 parameters
    params_tuple = ADMParams.getADMparams(self._R, self._T_base, self._T_ad)

    # Get substrate-dependent parameters
    substrate_params = self._get_substrate_dependent_params()

    # Unpack frequently used parameters
    (
        N_xc,
        N_I,
        N_aa,
        C_xc,
        C_sI,
        C_ch,
        C_pr,
        C_li,
        C_xI,
        C_su,
        C_aa,
        f_fa_li,
        C_fa,
        f_h2_su,
        f_bu_su,
        f_pro_su,
        f_ac_su,
        N_bac,
        C_bu,
        C_pro,
        C_ac,
        C_bac,
        Y_su,
        f_h2_aa,
        f_va_aa,
        f_bu_aa,
        f_pro_aa,
        f_ac_aa,
        C_va,
        Y_aa,
        Y_fa,
        Y_c4,
        Y_pro,
        C_ch4,
        Y_ac,
        Y_h2,
    ) = params_tuple[:36]

    # Get kinetic parameters (starting from index 36)
    kinetic_params = {
        "K_S_IN": params_tuple[36],
        "k_m_su": params_tuple[37],
        "K_S_su": params_tuple[38],
        "k_m_aa": params_tuple[41],
        "K_S_aa": params_tuple[42],
        "k_m_fa": params_tuple[43],
        "K_S_fa": params_tuple[44],
        "K_I_h2_fa": params_tuple[45],
        "k_m_c4": params_tuple[46],
        "K_S_c4": params_tuple[47],
        "K_I_h2_c4": params_tuple[48],
        "k_m_pro": params_tuple[49],
        "K_S_pro": params_tuple[50],
        "K_I_h2_pro": params_tuple[51],
        "k_m_ac": params_tuple[52],
        "K_S_ac": params_tuple[53],
        "K_I_nh3": params_tuple[54],
        "k_m_h2": params_tuple[57],
        "K_S_h2": params_tuple[58],
        "k_dec_X_su": params_tuple[61],
        "k_dec_X_aa": params_tuple[62],
        "k_dec_X_fa": params_tuple[63],
        "k_dec_X_c4": params_tuple[64],
        "k_dec_X_pro": params_tuple[65],
        "k_dec_X_ac": params_tuple[66],
        "k_dec_X_h2": params_tuple[67],
    }

    # Get acid-base and gas parameters
    acid_base_params = {
        "K_w": params_tuple[68],
        "K_a_va": params_tuple[69],
        "K_a_bu": params_tuple[70],
        "K_a_pro": params_tuple[71],
        "K_a_ac": params_tuple[72],
        "K_a_co2": params_tuple[73],
        "K_a_IN": params_tuple[74],
        "k_A_B_va": params_tuple[75],
        "k_A_B_bu": params_tuple[76],
        "k_A_B_pro": params_tuple[77],
        "k_A_B_ac": params_tuple[78],
        "k_A_B_co2": params_tuple[79],
        "k_A_B_IN": params_tuple[80],
    }

    gas_params = {
        "k_p": params_tuple[82],
        "k_L_a": params_tuple[83],
        "K_H_co2": params_tuple[84],
        "K_H_ch4": params_tuple[85],
        "K_H_h2": params_tuple[86],
    }

    # Get pH and inhibition parameters
    K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2 = ADMParams.getADMinhibitionparams()

    # Calculate pH and H+ concentration
    pH = ADMstate.calcPHOfADMstate(state_zero)
    S_H_ion = 10 ** (-pH)

    # Unpack state variables
    S_su, S_aa, S_fa, S_va, S_bu, S_pro, S_ac, S_h2 = state_zero[0:8]
    S_ch4, S_co2, S_nh4_ion, S_I = state_zero[8:12]
    X_xc, X_ch, X_pr, X_li = state_zero[12:16]
    X_su, X_aa, X_fa, X_c4, X_pro, X_ac, X_h2, X_I, X_p = state_zero[16:25]
    S_cation, S_anion = state_zero[25:27]
    S_va_ion, S_bu_ion, S_pro_ion, S_ac_ion = state_zero[27:31]
    S_hco3_ion, S_nh3 = state_zero[31:33]
    p_gas_h2, p_gas_ch4, p_gas_co2, pTOTAL = state_zero[33:37]

    # Get influent concentrations
    q_ad = np.sum(self._Q) if self._Q is not None else 0.0
    state_in = self._state_input if self._state_input is not None else [0.0] * 34

    # Calculate inhibition factors
    bio = BiochemicalProcesses()
    inhibitions = bio.calculate_inhibition_factors(
        S_H_ion,
        S_h2,
        S_nh4_ion,
        S_nh3,
        K_pH_aa,
        nn_aa,
        K_pH_ac,
        n_ac,
        K_pH_h2,
        n_h2,
        kinetic_params["K_S_IN"],
        kinetic_params["K_I_h2_fa"],
        kinetic_params["K_I_h2_c4"],
        kinetic_params["K_I_h2_pro"],
        kinetic_params["K_I_nh3"],
    )

    # Calculate biochemical process rates
    process_rates = bio.calculate_process_rates(state_zero, inhibitions, kinetic_params, substrate_params)
    Rho_1, Rho_2, Rho_3, Rho_4, Rho_5, Rho_6, Rho_7, Rho_8, Rho_9 = process_rates[:9]
    Rho_10, Rho_11, Rho_12, Rho_13, Rho_14, Rho_15, Rho_16 = process_rates[9:16]
    Rho_17, Rho_18, Rho_19 = process_rates[16:19]

    # Calculate acid-base rates
    acid_base_rates = bio.calculate_acid_base_rates(state_zero, acid_base_params)
    Rho_A_4, Rho_A_5, Rho_A_6, Rho_A_7, Rho_A_10, Rho_A_11 = acid_base_rates

    # Calculate gas transfer rates
    gas_rates = bio.calculate_gas_transfer_rates(state_zero, gas_params, self._RT, self.V_liq, self._V_gas, self._p_ext)
    Rho_T_8, Rho_T_9, Rho_T_10, Rho_T_11 = gas_rates

    # Extract substrate-dependent fractions
    f_ch_xc, f_pr_xc, f_li_xc, f_xI_xc, f_sI_xc, f_xp_xc, _, _, _, _, _, _, _, _ = substrate_params.values()

    # Calculate biomass decay fractions
    f_p = 0.08
    f_ch_xb = f_ch_xc / (f_ch_xc + f_pr_xc + f_li_xc) * (1 - f_p)
    f_pr_xb = f_pr_xc / (f_ch_xc + f_pr_xc + f_li_xc) * (1 - f_p)
    f_li_xb = f_li_xc / (f_ch_xc + f_pr_xc + f_li_xc) * (1 - f_p)

    # Differential equations for soluble components (1-12)
    diff_S_su = q_ad / self.V_liq * (state_in[0] - S_su) + Rho_2 + (1 - f_fa_li) * Rho_4 - Rho_5

    diff_S_aa = q_ad / self.V_liq * (state_in[1] - S_aa) + Rho_3 - Rho_6

    diff_S_fa = q_ad / self.V_liq * (state_in[2] - S_fa) + f_fa_li * Rho_4 - Rho_7

    diff_S_va = q_ad / self.V_liq * (state_in[3] - S_va) + (1 - Y_aa) * f_va_aa * Rho_6 - Rho_8

    diff_S_bu = (
        q_ad / self.V_liq * (state_in[4] - S_bu) + (1 - Y_su) * f_bu_su * Rho_5 + (1 - Y_aa) * f_bu_aa * Rho_6 - Rho_9
    )

    diff_S_pro = (
        q_ad / self.V_liq * (state_in[5] - S_pro)
        + (1 - Y_su) * f_pro_su * Rho_5
        + (1 - Y_aa) * f_pro_aa * Rho_6
        + (1 - Y_c4) * 0.54 * Rho_8
        - Rho_10
    )

    diff_S_ac = (
        q_ad / self.V_liq * (state_in[6] - S_ac)
        + (1 - Y_su) * f_ac_su * Rho_5
        + (1 - Y_aa) * f_ac_aa * Rho_6
        + (1 - Y_fa) * 0.7 * Rho_7
        + (1 - Y_c4) * 0.31 * Rho_8
        + (1 - Y_c4) * 0.8 * Rho_9
        + (1 - Y_pro) * 0.57 * Rho_10
        - Rho_11
    )

    diff_S_h2 = (
        q_ad / self.V_liq * (state_in[7] - S_h2)
        + (1 - Y_su) * f_h2_su * Rho_5
        + (1 - Y_aa) * f_h2_aa * Rho_6
        + (1 - Y_fa) * 0.3 * Rho_7
        + (1 - Y_c4) * 0.15 * Rho_8
        + (1 - Y_c4) * 0.2 * Rho_9
        + (1 - Y_pro) * 0.43 * Rho_10
        - Rho_12
        - self._V_gas / self.V_liq * Rho_T_8
    )

    diff_S_ch4 = (
        q_ad / self.V_liq * (state_in[8] - S_ch4)
        + (1 - Y_ac) * Rho_11
        + (1 - Y_h2) * Rho_12
        - self._V_gas / self.V_liq * Rho_T_9
    )

    # CO2 balance with carbon stoichiometry
    s_1 = -C_xc + f_sI_xc * C_sI + f_ch_xc * C_ch + f_pr_xc * C_pr + f_li_xc * C_li + f_xI_xc * C_xI
    s_2 = -C_ch + C_su
    s_3 = -C_pr + C_aa
    s_4 = -C_li + (1 - f_fa_li) * C_su + f_fa_li * C_fa
    s_5 = -C_su + (1 - Y_su) * (f_bu_su * C_bu + f_pro_su * C_pro + f_ac_su * C_ac) + Y_su * C_bac
    s_6 = -C_aa + (1 - Y_aa) * (f_va_aa * C_va + f_bu_aa * C_bu + f_pro_aa * C_pro + f_ac_aa * C_ac) + Y_aa * C_bac
    s_7 = -C_fa + (1 - Y_fa) * 0.7 * C_ac + Y_fa * C_bac
    s_8 = -C_va + (1 - Y_c4) * 0.54 * C_pro + (1 - Y_c4) * 0.31 * C_ac + Y_c4 * C_bac
    s_9 = -C_bu + (1 - Y_c4) * 0.8 * C_ac + Y_c4 * C_bac
    s_10 = -C_pro + (1 - Y_pro) * 0.57 * C_ac + Y_pro * C_bac
    s_11 = -C_ac + (1 - Y_ac) * C_ch4 + Y_ac * C_bac
    s_12 = (1 - Y_h2) * C_ch4 + Y_h2 * C_bac
    s_13 = -C_bac + C_xc

    Sigma = (
        s_1 * Rho_1
        + s_2 * Rho_2
        + s_3 * Rho_3
        + s_4 * Rho_4
        + s_5 * Rho_5
        + s_6 * Rho_6
        + s_7 * Rho_7
        + s_8 * Rho_8
        + s_9 * Rho_9
        + s_10 * Rho_10
        + s_11 * Rho_11
        + s_12 * Rho_12
        + s_13 * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
    )

    diff_S_co2 = q_ad / self.V_liq * (state_in[9] - S_co2) - Sigma - self._V_gas / self.V_liq * Rho_T_10 + Rho_A_10

    # Nitrogen balance
    diff_S_nh4_ion = (
        q_ad / self.V_liq * (state_in[10] - S_nh4_ion)
        - Y_su * N_bac * Rho_5
        + (N_aa - Y_aa * N_bac) * Rho_6
        - Y_fa * N_bac * Rho_7
        - Y_c4 * N_bac * Rho_8
        - Y_c4 * N_bac * Rho_9
        - Y_pro * N_bac * Rho_10
        - Y_ac * N_bac * Rho_11
        - Y_h2 * N_bac * Rho_12
        + (N_bac - N_xc) * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
        + Rho_A_11
    )

    diff_S_I = q_ad / self.V_liq * (state_in[11] - S_I) + f_sI_xc * Rho_1

    # Differential equations for particulate components (13-24)
    diff_X_xc = q_ad / self.V_liq * (state_in[12] - X_xc) - Rho_1

    diff_X_ch = (
        q_ad / self.V_liq * (state_in[13] - X_ch)
        + f_ch_xc * Rho_1
        - Rho_2
        + f_ch_xb * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
    )

    diff_X_pr = (
        q_ad / self.V_liq * (state_in[14] - X_pr)
        + f_pr_xc * Rho_1
        - Rho_3
        + f_pr_xb * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
    )

    diff_X_li = (
        q_ad / self.V_liq * (state_in[15] - X_li)
        + f_li_xc * Rho_1
        - Rho_4
        + f_li_xb * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
    )

    diff_X_su = q_ad / self.V_liq * (state_in[16] - X_su) + Y_su * Rho_5 - Rho_13
    diff_X_aa = q_ad / self.V_liq * (state_in[17] - X_aa) + Y_aa * Rho_6 - Rho_14
    diff_X_fa = q_ad / self.V_liq * (state_in[18] - X_fa) + Y_fa * Rho_7 - Rho_15
    diff_X_c4 = q_ad / self.V_liq * (state_in[19] - X_c4) + Y_c4 * Rho_8 + Y_c4 * Rho_9 - Rho_16
    diff_X_pro = q_ad / self.V_liq * (state_in[20] - X_pro) + Y_pro * Rho_10 - Rho_17
    diff_X_ac = q_ad / self.V_liq * (state_in[21] - X_ac) + Y_ac * Rho_11 - Rho_18
    diff_X_h2 = q_ad / self.V_liq * (state_in[22] - X_h2) + Y_h2 * Rho_12 - Rho_19
    diff_X_I = q_ad / self.V_liq * (state_in[23] - X_I) + f_xI_xc * Rho_1

    diff_X_p = (
        q_ad / self.V_liq * (state_in[24] - X_p)
        - f_xp_xc * Rho_1
        + f_p * (Rho_13 + Rho_14 + Rho_15 + Rho_16 + Rho_17 + Rho_18 + Rho_19)
    )

    # Differential equations for ions (25-32)
    diff_S_cation = q_ad / self.V_liq * (state_in[25] - S_cation)
    diff_S_anion = q_ad / self.V_liq * (state_in[26] - S_anion)
    diff_S_va_ion = q_ad / self.V_liq * (state_in[27] - S_va_ion) - Rho_A_4
    diff_S_bu_ion = q_ad / self.V_liq * (state_in[28] - S_bu_ion) - Rho_A_5
    diff_S_pro_ion = q_ad / self.V_liq * (state_in[29] - S_pro_ion) - Rho_A_6
    diff_S_ac_ion = q_ad / self.V_liq * (state_in[30] - S_ac_ion) - Rho_A_7
    diff_S_hco3_ion = q_ad / self.V_liq * (state_in[31] - S_hco3_ion) - Rho_A_10
    diff_S_nh3 = q_ad / self.V_liq * (state_in[32] - S_nh3) - Rho_A_11

    # Differential equations for gas phase (33-36)
    diff_p_gas_h2 = Rho_T_8 * self._RT / 16 - p_gas_h2 / pTOTAL * Rho_T_11
    diff_p_gas_ch4 = Rho_T_9 * self._RT / 64 - p_gas_ch4 / pTOTAL * Rho_T_11
    diff_p_gas_co2 = Rho_T_10 * self._RT - p_gas_co2 / pTOTAL * Rho_T_11
    diff_pTOTAL = self._RT / 16 * Rho_T_8 + self._RT / 64 * Rho_T_9 + self._RT * Rho_T_10 - 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_ion,
        diff_S_I,
        diff_X_xc,
        diff_X_ch,
        diff_X_pr,
        diff_X_li,
        diff_X_su,
        diff_X_aa,
        diff_X_fa,
        diff_X_c4,
        diff_X_pro,
        diff_X_ac,
        diff_X_h2,
        diff_X_I,
        diff_X_p,
        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_ion,
        diff_S_nh3,
        diff_p_gas_h2,
        diff_p_gas_ch4,
        diff_p_gas_co2,
        diff_pTOTAL,
    )

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

Initialize ADM1 model.

Parameters:

Name	Type	Description	Default
feedstock	Feedstock	Feedstock object containing substrate information. E.g. used to calculate ADM1 input stream.	required
V_liq	float	Liquid volume [m³]	1977.0
V_gas	float	Gas volume [m³]	304.0
T_ad	float	Operating temperature [K] (default: 308.15 = 35°C)	308.15

Source code in pyadm1/core/adm1.py

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

    Args:
        feedstock: Feedstock object containing substrate information. E.g. used to calculate ADM1 input stream.
        V_liq: Liquid volume [m³]
        V_gas: Gas volume [m³]
        T_ad: Operating temperature [K] (default: 308.15 = 35°C)
    """
    # Physical parameters
    self.V_liq = V_liq  # liquid volume of digester
    self._V_gas = V_gas  # gas volume of digester
    self._V_ad = self.V_liq + self._V_gas  # total volume of digester: liquid + gas volume
    self._T_ad = T_ad  # temperature inside the digester

    # Constants
    self._R = 0.08314  # 0.083145  Gas constant [bar·M^-1·K^-1]
    self._T_base = 298.15  # Base temperature [K] 25°C
    self._p_atm = 1.013  # Atmospheric pressure [bar]

    # Calculated parameters
    self._RT = self._R * self._T_ad  # R * T_ad
    # external pressures
    self._p_ext = self._p_atm - 0.0084147 * np.exp(0.054 * (self._T_ad - 273.15))

    # Feedstock and state
    self._feedstock = feedstock
    # vector of volumetric flow rates of the substrates. Length must be equal to the number of
    # substrates defined in xml
    self._Q: Optional[List[float]] = None
    self._state_input: Optional[List[float]] = None  # contains ADM1 input stream as a 34dim vector

    # Result tracking lists
    # Result tracking lists
    self._Q_GAS: List[float] = []  # produced biogas over all simulations [m³/d]
    self._Q_CH4: List[float] = []  # produced methane over all simulations [m³/d]
    self._Q_CO2: List[float] = []  # produced CO2 over all simulations [m³/d]
    self._P_GAS: List[float] = []  # gas pressures over all simulations [bar]
    self._pH_l: List[float] = []  # pH values over all simulations [-]
    self._FOSTAC: List[float] = []  # ratio of VFA over TA over all simulations [-]
    self._AcvsPro: List[float] = []  # ratio of acetic over propionic acid over all simulations [-]
    self._VFA: List[float] = []  # VFA concentrations over all simulations [g/L]
    self._TAC: List[float] = []  # TA concentrations over all simulations [g CaCO3/L]

calc_gas(pi_Sh2, pi_Sch4, pi_Sco2, pTOTAL)

Calculate biogas production rates from partial pressures.

Uses the ideal gas law and Henry's constants to calculate gas flow rates from the gas phase partial pressures.

Parameters:

Name	Type	Description	Default
pi_Sh2	float	Hydrogen partial pressure [bar]	required
pi_Sch4	float	Methane partial pressure [bar]	required
pi_Sco2	float	CO2 partial pressure [bar]	required
pTOTAL	float	Total gas pressure [bar]	required

Returns:

Type	Description
Tuple[float, float, float, float]	Tuple containing: - q_gas: Total biogas flow rate [m³/d] - q_ch4: Methane flow rate [m³/d] - q_co2: CO2 flow rate [m³/d] - p_gas: Total gas partial pressure (excl. H2O) [bar]

Example

q_gas, q_ch4, q_co2, p_gas = adm1.calc_gas(5e-6, 0.55, 0.42, 0.98) print(f"Biogas: {q_gas:.1f} m³/d, Methane: {q_ch4:.1f} m³/d")

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]:
    """
    Calculate biogas production rates from partial pressures.

    Uses the ideal gas law and Henry's constants to calculate gas flow rates
    from the gas phase partial pressures.

    Args:
        pi_Sh2: Hydrogen partial pressure [bar]
        pi_Sch4: Methane partial pressure [bar]
        pi_Sco2: CO2 partial pressure [bar]
        pTOTAL: Total gas pressure [bar]

    Returns:
        Tuple containing:
            - q_gas: Total biogas flow rate [m³/d]
            - q_ch4: Methane flow rate [m³/d]
            - q_co2: CO2 flow rate [m³/d]
            - p_gas: Total gas partial pressure (excl. H2O) [bar]

    Example:
        >>> q_gas, q_ch4, q_co2, p_gas = adm1.calc_gas(5e-6, 0.55, 0.42, 0.98)
        >>> print(f"Biogas: {q_gas:.1f} m³/d, Methane: {q_ch4:.1f} m³/d")
    """
    _, k_p, _, _, _, _ = ADMParams.getADMgasparams(self._R, self._T_base, self._T_ad)

    # Ideal gas law constant for conversion
    NQ = 44.643

    # Total biogas flow from pressure difference
    q_gas = k_p * (pTOTAL - self._p_ext) / (self._RT / 1000 * NQ) * self.V_liq

    # Total gas partial pressure (excluding water vapor)
    p_gas = pi_Sh2 + pi_Sch4 + pi_Sco2

    # Ensure non-negative gas flows
    if isinstance(q_gas, np.ndarray):
        q_gas = np.maximum(q_gas, 0.0)
    else:
        q_gas = max(q_gas, 0.0)

    # Calculate component flows from partial pressures
    if p_gas > 0:
        q_ch4 = q_gas * (pi_Sch4 / p_gas)
        q_co2 = q_gas * (pi_Sco2 / p_gas)
    else:
        q_ch4 = 0.0
        q_co2 = 0.0

    # Ensure non-negative
    if isinstance(q_ch4, np.ndarray):
        q_ch4 = np.maximum(q_ch4, 0.0)
        q_co2 = np.maximum(q_co2, 0.0)
    else:
        q_ch4 = max(q_ch4, 0.0)
        q_co2 = max(q_co2, 0.0)

    return q_gas, q_ch4, q_co2, p_gas

clear_calibration_parameters()

Clear all calibration parameters and revert to substrate-dependent calculations.

Example

adm1.clear_calibration_parameters()

Source code in pyadm1/core/adm1.py

def clear_calibration_parameters(self) -> None:
    """
    Clear all calibration parameters and revert to substrate-dependent calculations.

    Example:
        >>> adm1.clear_calibration_parameters()
    """
    if hasattr(self, "_calibration_params"):
        self._calibration_params = {}

create_influent(Q, i)

Create ADM1 input stream from volumetric flow rates.

Calculates the ADM1 influent state by mixing substrate streams according to their volumetric flow rates. The resulting influent composition is stored internally for use in ODE calculations.

Parameters:

Name	Type	Description	Default
Q	List[float]	Volumetric flow rates for each substrate [m³/d] Length must equal number of substrates in feedstock	required
i	int	Time step index for accessing influent dataframe	required

Example

adm1.create_influent([15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 0)

Source code in pyadm1/core/adm1.py

def create_influent(self, Q: List[float], i: int) -> None:
    """
    Create ADM1 input stream from volumetric flow rates.

    Calculates the ADM1 influent state by mixing substrate streams according
    to their volumetric flow rates. The resulting influent composition is
    stored internally for use in ODE calculations.

    Args:
        Q: Volumetric flow rates for each substrate [m³/d]
            Length must equal number of substrates in feedstock
        i: Time step index for accessing influent dataframe

    Example:
        >>> adm1.create_influent([15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 0)
    """
    self._Q = Q
    influent_state = self._feedstock.get_influent_dataframe(Q)
    self._set_influent(influent_state, i)

get_calibration_parameters()

Get currently set calibration parameters.

Returns:

Name	Type	Description
dict	dict	Current calibration parameters as {param_name: value}.

Example

params = adm1.get_calibration_parameters() print(params)

Source code in pyadm1/core/adm1.py

def get_calibration_parameters(self) -> dict:
    """
    Get currently set calibration parameters.

    Returns:
        dict: Current calibration parameters as {param_name: value}.

    Example:
        >>> params = adm1.get_calibration_parameters()
        >>> print(params)
        {'k_dis': 0.55, 'Y_su': 0.105}
    """
    if hasattr(self, "_calibration_params"):
        return self._calibration_params.copy()
    return {}

print_params_at_current_state(state_ADM1xp)

Calculate and store process parameters from current state.

Computes key process indicators including pH, VFA, TAC, and gas production rates. Also stores values in tracking lists.

Parameters:

Name	Type	Description	Default
state_ADM1xp	List[float]	Current ADM1 state vector (37 elements)	required

Example

adm1.print_params_at_current_state(state_vector)

Source code in pyadm1/core/adm1.py

def print_params_at_current_state(self, state_ADM1xp: List[float]) -> None:
    """
    Calculate and store process parameters from current state.

    Computes key process indicators including pH, VFA, TAC,
    and gas production rates. Also stores values in tracking lists.

    Args:
        state_ADM1xp: Current ADM1 state vector (37 elements)

    Example:
        >>> adm1.print_params_at_current_state(state_vector)
    """
    # Calculate process indicators using DLL
    self._pH_l.append(np.round(ADMstate.calcPHOfADMstate(state_ADM1xp), 1))
    self._FOSTAC.append(np.round(ADMstate.calcFOSTACOfADMstate(state_ADM1xp).Value, 2))
    self._AcvsPro.append(np.round(ADMstate.calcAcetic_vs_PropionicOfADMstate(state_ADM1xp).Value, 1))
    self._VFA.append(np.round(ADMstate.calcVFAOfADMstate(state_ADM1xp, "gHAceq/l").Value, 2))
    self._TAC.append(np.round(ADMstate.calcTACOfADMstate(state_ADM1xp, "gCaCO3eq/l").Value, 1))

    # Ensure at least 2 values in lists (because the last three values go to the controller)
    # I am assuming here that we start from a steady state
    if len(self._pH_l) < 2:
        self._pH_l.append(self._pH_l[-1])
        self._FOSTAC.append(self._FOSTAC[-1])
        self._AcvsPro.append(self._AcvsPro[-1])
        self._VFA.append(self._VFA[-1])
        self._TAC.append(self._TAC[-1])

    # Calculate and store gas production
    q_gas, q_ch4, q_co2, p_gas = self.calc_gas(state_ADM1xp[33], state_ADM1xp[34], state_ADM1xp[35], state_ADM1xp[36])

    self._Q_GAS.append(q_gas)
    self._Q_CH4.append(q_ch4)
    self._Q_CO2.append(q_co2)
    self._P_GAS.append(p_gas)

    # Ensure at least 2 values, because the last three values go to controller.
    # I am assuming here that we start from a steady state
    if len(self._Q_GAS) < 2:
        for _ in range(2):  # to have at least 4 values in the lists
            self._Q_GAS.append(q_gas)
            self._Q_CH4.append(q_ch4)
            self._Q_CO2.append(q_co2)
            self._P_GAS.append(p_gas)

save_final_state_in_csv(simulate_results, filename='digester_final.csv')

Save final ADM1 state vector to CSV file.

Exports only the last state from simulation results, which can be used as initial state for subsequent simulations.

Parameters:

Name	Type	Description	Default
simulate_results	List[List[float]]	List of ADM1 state vectors from simulation	required
filename	str	Output CSV filename	'digester_final.csv'

Example

results = [[0.01]37, [0.02]37, [0.03]*37] adm1.save_final_state_in_csv(results, 'final_state.csv')

Source code in pyadm1/core/adm1.py

def save_final_state_in_csv(self, simulate_results: List[List[float]], filename: str = "digester_final.csv") -> None:
    """
    Save final ADM1 state vector to CSV file.

    Exports only the last state from simulation results, which can be used
    as initial state for subsequent simulations.

    Args:
        simulate_results: List of ADM1 state vectors from simulation
        filename: Output CSV filename

    Example:
        >>> results = [[0.01]*37, [0.02]*37, [0.03]*37]
        >>> adm1.save_final_state_in_csv(results, 'final_state.csv')
    """
    columns = [
        *self._feedstock.header()[:-1],
        "pi_Sh2",
        "pi_Sch4",
        "pi_Sco2",
        "pTOTAL",
    ]

    simulate_results_df = pd.DataFrame.from_records(simulate_results)
    simulate_results_df.columns = columns

    # Keep only the last row
    last_simulated_result = simulate_results_df[-1:]
    last_simulated_result.to_csv(filename, index=False)

set_calibration_parameters(parameters)

Set calibration parameters that override substrate-dependent calculations.

Parameters:

Name	Type	Description	Default
parameters	dict	Parameter values as {param_name: value}.	required

Example

adm1.set_calibration_parameters({ ... 'k_dis': 0.55, ... 'k_hyd_ch': 11.0, ... 'Y_su': 0.105 ... })

Source code in pyadm1/core/adm1.py

def set_calibration_parameters(self, parameters: dict) -> None:
    """
    Set calibration parameters that override substrate-dependent calculations.

    Args:
        parameters: Parameter values as {param_name: value}.

    Example:
        >>> adm1.set_calibration_parameters({
        ...     'k_dis': 0.55,
        ...     'k_hyd_ch': 11.0,
        ...     'Y_su': 0.105
        ... })
    """
    if not hasattr(self, "_calibration_params"):
        self._calibration_params = {}

    self._calibration_params.update(parameters)

pyadm1.core.adm_params.ADMParams

Static class containing ADM1 model parameters.

Source code in pyadm1/core/adm_params.py

class ADMParams:
    """Static class containing ADM1 model parameters."""

    # *** CONSTRUCTORS ***

    # *** PUBLIC SET methods ***

    # *** PUBLIC GET methods ***

    # *** PUBLIC methods ***

    # *** PUBLIC STATIC/CLASS GET methods ***

    @staticmethod
    def getADMparams(R: float, T_base: float, T_ad: float) -> Tuple[float, ...]:
        """
        Get all ADM1 stoichiometric and kinetic parameters.

        Parameters
        ----------
        R : float
            Gas constant [bar·M^-1·K^-1]
        T_base : float
            Base temperature [K]
        T_ad : float
            Digester temperature [K]

        Returns
        -------
        Tuple[float, ...]
            All ADM1 parameters (87 values)
        """
        # parameter definition from the Rosen et al (2006) BSM2 report bmadm1_report
        # Stoichiometric parameter
        # they are substrate dependent and therefore calculated from the current substrate mix, directly in PyADM1.py
        # f_sI_xc =  0.1      # OK
        # f_xI_xc =  0.2      # in C# split into fXI_XC and fXP_XC
        # f_ch_xc =  0.2      # OK
        # f_pr_xc =  0.2      # OK
        # f_li_xc =  0.3      # OK

        N_xc = 0.0376 / 14  # OK
        N_I = 0.06 / 14  # kmole N.kg^-1COD  # OK
        N_aa = 0.007  # kmole N.kg^-1COD         OK
        C_xc = 0.03  # kmole C.kg^-1COD
        C_sI = 0.03  # kmole C.kg^-1COD          OK
        C_ch = 0.0313  # kmole C.kg^-1COD        C_Xch OK
        C_pr = 0.03  # kmole C.kg^-1COD          C_Xpr OK
        C_li = 0.022  # kmole C.kg^-1COD         C_Xli OK
        C_xI = 0.03  # kmole C.kg^-1COD          OK
        C_su = 0.0313  # kmole C.kg^-1COD        OK
        C_aa = 0.03  # kmole C.kg^-1COD          OK
        f_fa_li = 0.95  # fFA_Xli OK
        C_fa = 0.0217  # kmole C.kg^-1COD        # C_Sfa OK

        f_h2_su, f_bu_su, f_pro_su, f_ac_su = ADMParams._getADMfsuparams()

        N_bac = 0.08 / 14  # kmole N.kg^-1COD        # N_XB OK
        C_bu = 0.025  # kmole C.kg^-1COD         C_Sbu OK
        C_pro = 0.0268  # kmole C.kg^-1COD       C_Spro OK
        C_ac = 0.0313  # kmole C.kg^-1COD        C_Sac OK
        C_bac = 0.0313  # kmole C.kg^-1COD       C_XB OK

        C_va = 0.024  # kmole C.kg^-1COD        # C_Sva OK
        C_ch4 = 0.0156  # kmole C.kg^-1COD      # C_Sch4

        Y_su, Y_aa, Y_fa, Y_c4, Y_pro, Y_ac, Y_h2 = ADMParams._getADMYparams()

        f_h2_aa, f_va_aa, f_bu_aa, f_pro_aa, f_ac_aa = ADMParams._getADMfaaparams()

        # they are substrate dependent and therefore calculated from the current substrate mix, directly in PyADM1.py
        # Biochemical parameter values from the Rosen et al (2006) BSM2 report
        # k_dis =  0.5 #d^-1              OK
        # k_hyd_ch =  10 #d^-1            OK
        # k_hyd_pr =  10 #d^-1            OK
        # k_hyd_li =  10 #d^-1            OK

        (
            K_S_IN,
            k_m_su,
            K_S_su,
            k_m_aa,
            K_S_aa,
            k_m_fa,
            K_S_fa,
            K_I_h2_fa,
            k_m_c4,
            K_S_c4,
            K_I_h2_c4,
            k_m_pro,
            K_S_pro,
            K_I_h2_pro,
            k_m_ac,
            K_S_ac,
            K_I_nh3,
            k_m_h2,
            K_S_h2,
        ) = ADMParams._getADMk_mK_Sparams()

        pH_LL_aa, pH_UL_aa, pH_LL_ac, pH_UL_ac, pH_LL_h2, pH_UL_h2 = ADMParams._getADMpHULLLparams()

        # decay rates
        k_dec_X_su = 0.02  # d^-1            # OK
        k_dec_X_aa = 0.02  # d^-1            # OK
        k_dec_X_fa = 0.02  # d^-1            # OK
        k_dec_X_c4 = 0.02  # d^-1            # OK
        k_dec_X_pro = 0.02  # d^-1           # OK
        k_dec_X_ac = 0.02  # d^-1            # OK
        k_dec_X_h2 = 0.02  # d^-1            # OK
        ## M is kmole m^-3

        K_w, K_a_va, K_a_bu, K_a_pro, K_a_ac, K_a_co2, K_a_IN = ADMParams._getADMKparams(R, T_base, T_ad)

        # acid-base kinetic parameters
        # those values all are 1e8 kmole/d in C# implementation
        k_A_B_va = 1e8  # 10 ** 10 #M^-1 * d^-1
        k_A_B_bu = 1e8  # 10 ** 10 #M^-1 * d^-1
        k_A_B_pro = 1e8  # 10 ** 10 #M^-1 * d^-1
        k_A_B_ac = 1e8  # 10 ** 10 #M^-1 * d^-1
        k_A_B_co2 = 1e8  # 10 ** 10 #M^-1 * d^-1
        k_A_B_IN = 1e8  # 10 ** 10 #M^-1 * d^-1

        p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2 = ADMParams.getADMgasparams(R, T_base, T_ad)

        return (
            N_xc,
            N_I,
            N_aa,
            C_xc,
            C_sI,
            C_ch,
            C_pr,
            C_li,
            C_xI,
            C_su,
            C_aa,
            f_fa_li,
            C_fa,
            f_h2_su,
            f_bu_su,
            f_pro_su,
            f_ac_su,
            N_bac,
            C_bu,
            C_pro,
            C_ac,
            C_bac,
            Y_su,
            f_h2_aa,
            f_va_aa,
            f_bu_aa,
            f_pro_aa,
            f_ac_aa,
            C_va,
            Y_aa,
            Y_fa,
            Y_c4,
            Y_pro,
            C_ch4,
            Y_ac,
            Y_h2,
            K_S_IN,
            k_m_su,
            K_S_su,
            pH_UL_aa,
            pH_LL_aa,
            k_m_aa,
            K_S_aa,
            k_m_fa,
            K_S_fa,
            K_I_h2_fa,
            k_m_c4,
            K_S_c4,
            K_I_h2_c4,
            k_m_pro,
            K_S_pro,
            K_I_h2_pro,
            k_m_ac,
            K_S_ac,
            K_I_nh3,
            pH_UL_ac,
            pH_LL_ac,
            k_m_h2,
            K_S_h2,
            pH_UL_h2,
            pH_LL_h2,
            k_dec_X_su,
            k_dec_X_aa,
            k_dec_X_fa,
            k_dec_X_c4,
            k_dec_X_pro,
            k_dec_X_ac,
            k_dec_X_h2,
            K_w,
            K_a_va,
            K_a_bu,
            K_a_pro,
            K_a_ac,
            K_a_co2,
            K_a_IN,
            k_A_B_va,
            k_A_B_bu,
            k_A_B_pro,
            k_A_B_ac,
            k_A_B_co2,
            k_A_B_IN,
            p_gas_h2o,
            k_p,
            k_L_a,
            K_H_co2,
            K_H_ch4,
            K_H_h2,
        )

    @staticmethod
    def getADMinhibitionparams() -> Tuple[float, float, float, float, float, float]:
        """
        Get pH inhibition parameters.

        Returns
        -------
        Tuple[float, float, float, float, float, float]
            K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2
        """
        pH_LL_aa, pH_UL_aa, pH_LL_ac, pH_UL_ac, pH_LL_h2, pH_UL_h2 = ADMParams._getADMpHULLLparams()

        # related to pH inhibition taken from BSM2 report, they are global variables to avoid repeating them in DAE part
        K_pH_aa = 10 ** (-1 * (pH_LL_aa + pH_UL_aa) / 2.0)  # OK
        # we need a difference between N_aa and n_aa to avoid typos and nn_aa refers to the n_aa in BSM2 report
        # in C# just set to 2, the calculation here is also 2, so OK
        nn_aa = 3.0 / (pH_UL_aa - pH_LL_aa)
        K_pH_ac = 10 ** (-1 * (pH_LL_ac + pH_UL_ac) / 2.0)  # OK
        n_ac = 3.0 / (pH_UL_ac - pH_LL_ac)  # OK
        K_pH_h2 = 10 ** (-1 * (pH_LL_h2 + pH_UL_h2) / 2.0)  # OK
        n_h2 = 3.0 / (pH_UL_h2 - pH_LL_h2)  # OK

        return K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2

    @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^-1·K^-1]
        T_base : float
            Base temperature [K]
        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))  # bar #0.0557
        # in C# this value is 10000, but in m^3/(m^3*d)
        # k_p = 5 * 10 ** 4  # m^3.d^-1.bar^-1 #only for BSM2 AD conditions, recalibrate for other AD cases #gas outlet friction
        k_p = 1 * 10**4
        k_L_a = 200.0  # d^-1 OK in C# implementation there is klAH2, klaCH4, klaCO2, but all are equal to 200
        K_H_co2 = 1 / (0.0271 * R * T_ad)
        K_H_ch4 = 1 / (0.00116 * R * T_ad)
        K_H_h2 = 1 / (7.38e-4 * R * T_ad)

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

    # *** PRIVATE methods ***

    @staticmethod
    def _getADMKparams(R: float, T_base: float, T_ad: float) -> Tuple[float, float, float, float, float, float, float]:
        """
        Return acid-base equilibrium coefficients.

        Parameters
        ----------
        R : float
            Gas constant [bar·M^-1·K^-1]
        T_base : float
            Base temperature [K]
        T_ad : float
            Digester temperature [K]

        Returns
        -------
        Tuple[float, float, float, float, float, float, float]
            K_w, K_a_va, K_a_bu, K_a_pro, K_a_ac, K_a_co2, K_a_IN
        """
        K_w = 10**-14.0 * np.exp((55900 / (100 * R)) * (1 / T_base - 1 / T_ad))  # M #2.08 * 10 ^ -14  OK

        K_a_va = 10**-4.86  # M  ADM1 value = 1.38 * 10 ^ -5      OK
        K_a_bu = 10**-4.82  # M #1.5 * 10 ^ -5                    OK
        K_a_pro = 10**-4.88  # M #1.32 * 10 ^ -5                  OK
        K_a_ac = 10**-4.76  # M #1.74 * 10 ^ -5                   OK

        K_a_co2 = 10**-6.35 * np.exp((7646 / (100 * R)) * (1 / T_base - 1 / T_ad))  # M #4.94 * 10 ^ -7 OK
        K_a_IN = 10**-9.25 * np.exp((51965 / (100 * R)) * (1 / T_base - 1 / T_ad))  # M #1.11 * 10 ^ -9 OK

        return K_w, K_a_va, K_a_bu, K_a_pro, K_a_ac, K_a_co2, K_a_IN

    @staticmethod
    def _getADMpHULLLparams() -> Tuple[float, float, float, float, float, float]:
        """
        Return upper and lower limits for inhibition by pH.

        Returns
        -------
        Tuple[float, float, float, float, float, float]
            pH_LL_aa, pH_UL_aa, pH_LL_ac, pH_UL_ac, pH_LL_h2, pH_UL_h2
        """
        pH_UL_aa = 5.5  # OK
        pH_LL_aa = 4  # OK
        pH_UL_ac = 7  # OK
        pH_LL_ac = 6  # OK
        pH_UL_h2 = 6  # OK
        pH_LL_h2 = 5  # OK

        return pH_LL_aa, pH_UL_aa, pH_LL_ac, pH_UL_ac, pH_LL_h2, pH_UL_h2

    @staticmethod
    def _getADMk_mK_Sparams() -> Tuple[float, ...]:
        """
        Return maximum uptake rates and half-saturation constants.

        Returns
        -------
        Tuple[float, ...]
            K_S_IN, k_m_su, K_S_su, k_m_aa, K_S_aa, k_m_fa, K_S_fa, K_I_h2_fa,
            k_m_c4, K_S_c4, K_I_h2_c4, k_m_pro, K_S_pro, K_I_h2_pro, k_m_ac,
            K_S_ac, K_I_nh3, k_m_h2, K_S_h2
        """
        K_S_IN = 10**-4  # M              OK
        k_m_su = 30  # d^-1                 OK
        K_S_su = 0.5  # kgCOD.m^-3          OK
        k_m_aa = 50  # d^-1                 OK
        K_S_aa = 0.3  ##kgCOD.m^-3          OK
        k_m_fa = 6  # d^-1                  OK
        K_S_fa = 0.4  # kgCOD.m^-3          OK
        K_I_h2_fa = 5 * 10**-6  # kgCOD.m^-3      OK
        k_m_c4 = 20  # d^-1                     OK
        K_S_c4 = 0.2  # kgCOD.m^-3 (SIMBA5 default)
        K_I_h2_c4 = 10**-5  # kgCOD.m^-3      OK
        k_m_pro = 13  # d^-1                    OK
        K_S_pro = 0.1  # kgCOD.m^-3             OK
        K_I_h2_pro = 3.5 * 10**-6  # kgCOD.m^-3       OK
        k_m_ac = 8  # kgCOD.m^-3                    OK
        K_S_ac = 0.15  # kgCOD.m^-3                  OK
        K_I_nh3 = 0.0018  # M (SIMBA5 default)
        k_m_h2 = 35  # d^-1             OK
        K_S_h2 = 7 * 10**-6  # kgCOD.m^-3         OK

        return (
            K_S_IN,
            k_m_su,
            K_S_su,
            k_m_aa,
            K_S_aa,
            k_m_fa,
            K_S_fa,
            K_I_h2_fa,
            k_m_c4,
            K_S_c4,
            K_I_h2_c4,
            k_m_pro,
            K_S_pro,
            K_I_h2_pro,
            k_m_ac,
            K_S_ac,
            K_I_nh3,
            k_m_h2,
            K_S_h2,
        )

    @staticmethod
    def _getADMYparams() -> Tuple[float, float, float, float, float, float, float]:
        """
        Return yield uptake parameters.

        Returns
        -------
        Tuple[float, float, float, float, float, float, float]
            Y_su, Y_aa, Y_fa, Y_c4, Y_pro, Y_ac, Y_h2
        """
        Y_su = 0.1  # OK
        Y_aa = 0.08  # OK
        Y_fa = 0.06  # OK
        Y_c4 = 0.06  # OK
        Y_pro = 0.04  # OK
        Y_ac = 0.05  # OK
        Y_h2 = 0.06  # OK

        return (Y_su, Y_aa, Y_fa, Y_c4, Y_pro, Y_ac, Y_h2)

    @staticmethod
    def _getADMfaaparams() -> Tuple[float, float, float, float, float]:
        """
        Return fraction parameters from amino acids.

        Returns
        -------
        Tuple[float, float, float, float, float]
            f_h2_aa, f_va_aa, f_bu_aa, f_pro_aa, f_ac_aa
        """
        f_h2_aa = 0.06  # OK
        f_va_aa = 0.23  # OK
        f_bu_aa = 0.26  # OK
        f_pro_aa = 0.05  # OK
        f_ac_aa = 0.40  # OK

        return f_h2_aa, f_va_aa, f_bu_aa, f_pro_aa, f_ac_aa

    @staticmethod
    def _getADMfsuparams() -> Tuple[float, float, float, float]:
        """
        Return fraction parameters from sugars.

        Returns
        -------
        Tuple[float, float, float, float]
            f_h2_su, f_bu_su, f_pro_su, f_ac_su
        """
        f_h2_su = 0.19  # OK
        f_bu_su = 0.13  # OK
        f_pro_su = 0.27  # OK
        f_ac_su = 0.41  # OK

        return f_h2_su, f_bu_su, f_pro_su, f_ac_su

Functions

getADMgasparams(R, T_base, T_ad) staticmethod

Get gas phase parameters including Henry constants.

Parameters

R : float Gas constant [bar·M^-1·K^-1] T_base : float Base temperature [K] 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^-1·K^-1]
    T_base : float
        Base temperature [K]
    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))  # bar #0.0557
    # in C# this value is 10000, but in m^3/(m^3*d)
    # k_p = 5 * 10 ** 4  # m^3.d^-1.bar^-1 #only for BSM2 AD conditions, recalibrate for other AD cases #gas outlet friction
    k_p = 1 * 10**4
    k_L_a = 200.0  # d^-1 OK in C# implementation there is klAH2, klaCH4, klaCO2, but all are equal to 200
    K_H_co2 = 1 / (0.0271 * R * T_ad)
    K_H_ch4 = 1 / (0.00116 * R * T_ad)
    K_H_h2 = 1 / (7.38e-4 * R * T_ad)

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

getADMinhibitionparams() staticmethod

Get pH inhibition parameters.

Returns

Tuple[float, float, float, float, float, float] K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2

Source code in pyadm1/core/adm_params.py

@staticmethod
def getADMinhibitionparams() -> Tuple[float, float, float, float, float, float]:
    """
    Get pH inhibition parameters.

    Returns
    -------
    Tuple[float, float, float, float, float, float]
        K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2
    """
    pH_LL_aa, pH_UL_aa, pH_LL_ac, pH_UL_ac, pH_LL_h2, pH_UL_h2 = ADMParams._getADMpHULLLparams()

    # related to pH inhibition taken from BSM2 report, they are global variables to avoid repeating them in DAE part
    K_pH_aa = 10 ** (-1 * (pH_LL_aa + pH_UL_aa) / 2.0)  # OK
    # we need a difference between N_aa and n_aa to avoid typos and nn_aa refers to the n_aa in BSM2 report
    # in C# just set to 2, the calculation here is also 2, so OK
    nn_aa = 3.0 / (pH_UL_aa - pH_LL_aa)
    K_pH_ac = 10 ** (-1 * (pH_LL_ac + pH_UL_ac) / 2.0)  # OK
    n_ac = 3.0 / (pH_UL_ac - pH_LL_ac)  # OK
    K_pH_h2 = 10 ** (-1 * (pH_LL_h2 + pH_UL_h2) / 2.0)  # OK
    n_h2 = 3.0 / (pH_UL_h2 - pH_LL_h2)  # OK

    return K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2

getADMparams(R, T_base, T_ad) staticmethod

Get all ADM1 stoichiometric and kinetic parameters.

Parameters

R : float Gas constant [bar·M^-1·K^-1] T_base : float Base temperature [K] T_ad : float Digester temperature [K]

Returns

Tuple[float, ...] All ADM1 parameters (87 values)

Source code in pyadm1/core/adm_params.py

@staticmethod
def getADMparams(R: float, T_base: float, T_ad: float) -> Tuple[float, ...]:
    """
    Get all ADM1 stoichiometric and kinetic parameters.

    Parameters
    ----------
    R : float
        Gas constant [bar·M^-1·K^-1]
    T_base : float
        Base temperature [K]
    T_ad : float
        Digester temperature [K]

    Returns
    -------
    Tuple[float, ...]
        All ADM1 parameters (87 values)
    """
    # parameter definition from the Rosen et al (2006) BSM2 report bmadm1_report
    # Stoichiometric parameter
    # they are substrate dependent and therefore calculated from the current substrate mix, directly in PyADM1.py
    # f_sI_xc =  0.1      # OK
    # f_xI_xc =  0.2      # in C# split into fXI_XC and fXP_XC
    # f_ch_xc =  0.2      # OK
    # f_pr_xc =  0.2      # OK
    # f_li_xc =  0.3      # OK

    N_xc = 0.0376 / 14  # OK
    N_I = 0.06 / 14  # kmole N.kg^-1COD  # OK
    N_aa = 0.007  # kmole N.kg^-1COD         OK
    C_xc = 0.03  # kmole C.kg^-1COD
    C_sI = 0.03  # kmole C.kg^-1COD          OK
    C_ch = 0.0313  # kmole C.kg^-1COD        C_Xch OK
    C_pr = 0.03  # kmole C.kg^-1COD          C_Xpr OK
    C_li = 0.022  # kmole C.kg^-1COD         C_Xli OK
    C_xI = 0.03  # kmole C.kg^-1COD          OK
    C_su = 0.0313  # kmole C.kg^-1COD        OK
    C_aa = 0.03  # kmole C.kg^-1COD          OK
    f_fa_li = 0.95  # fFA_Xli OK
    C_fa = 0.0217  # kmole C.kg^-1COD        # C_Sfa OK

    f_h2_su, f_bu_su, f_pro_su, f_ac_su = ADMParams._getADMfsuparams()

    N_bac = 0.08 / 14  # kmole N.kg^-1COD        # N_XB OK
    C_bu = 0.025  # kmole C.kg^-1COD         C_Sbu OK
    C_pro = 0.0268  # kmole C.kg^-1COD       C_Spro OK
    C_ac = 0.0313  # kmole C.kg^-1COD        C_Sac OK
    C_bac = 0.0313  # kmole C.kg^-1COD       C_XB OK

    C_va = 0.024  # kmole C.kg^-1COD        # C_Sva OK
    C_ch4 = 0.0156  # kmole C.kg^-1COD      # C_Sch4

    Y_su, Y_aa, Y_fa, Y_c4, Y_pro, Y_ac, Y_h2 = ADMParams._getADMYparams()

    f_h2_aa, f_va_aa, f_bu_aa, f_pro_aa, f_ac_aa = ADMParams._getADMfaaparams()

    # they are substrate dependent and therefore calculated from the current substrate mix, directly in PyADM1.py
    # Biochemical parameter values from the Rosen et al (2006) BSM2 report
    # k_dis =  0.5 #d^-1              OK
    # k_hyd_ch =  10 #d^-1            OK
    # k_hyd_pr =  10 #d^-1            OK
    # k_hyd_li =  10 #d^-1            OK

    (
        K_S_IN,
        k_m_su,
        K_S_su,
        k_m_aa,
        K_S_aa,
        k_m_fa,
        K_S_fa,
        K_I_h2_fa,
        k_m_c4,
        K_S_c4,
        K_I_h2_c4,
        k_m_pro,
        K_S_pro,
        K_I_h2_pro,
        k_m_ac,
        K_S_ac,
        K_I_nh3,
        k_m_h2,
        K_S_h2,
    ) = ADMParams._getADMk_mK_Sparams()

    pH_LL_aa, pH_UL_aa, pH_LL_ac, pH_UL_ac, pH_LL_h2, pH_UL_h2 = ADMParams._getADMpHULLLparams()

    # decay rates
    k_dec_X_su = 0.02  # d^-1            # OK
    k_dec_X_aa = 0.02  # d^-1            # OK
    k_dec_X_fa = 0.02  # d^-1            # OK
    k_dec_X_c4 = 0.02  # d^-1            # OK
    k_dec_X_pro = 0.02  # d^-1           # OK
    k_dec_X_ac = 0.02  # d^-1            # OK
    k_dec_X_h2 = 0.02  # d^-1            # OK
    ## M is kmole m^-3

    K_w, K_a_va, K_a_bu, K_a_pro, K_a_ac, K_a_co2, K_a_IN = ADMParams._getADMKparams(R, T_base, T_ad)

    # acid-base kinetic parameters
    # those values all are 1e8 kmole/d in C# implementation
    k_A_B_va = 1e8  # 10 ** 10 #M^-1 * d^-1
    k_A_B_bu = 1e8  # 10 ** 10 #M^-1 * d^-1
    k_A_B_pro = 1e8  # 10 ** 10 #M^-1 * d^-1
    k_A_B_ac = 1e8  # 10 ** 10 #M^-1 * d^-1
    k_A_B_co2 = 1e8  # 10 ** 10 #M^-1 * d^-1
    k_A_B_IN = 1e8  # 10 ** 10 #M^-1 * d^-1

    p_gas_h2o, k_p, k_L_a, K_H_co2, K_H_ch4, K_H_h2 = ADMParams.getADMgasparams(R, T_base, T_ad)

    return (
        N_xc,
        N_I,
        N_aa,
        C_xc,
        C_sI,
        C_ch,
        C_pr,
        C_li,
        C_xI,
        C_su,
        C_aa,
        f_fa_li,
        C_fa,
        f_h2_su,
        f_bu_su,
        f_pro_su,
        f_ac_su,
        N_bac,
        C_bu,
        C_pro,
        C_ac,
        C_bac,
        Y_su,
        f_h2_aa,
        f_va_aa,
        f_bu_aa,
        f_pro_aa,
        f_ac_aa,
        C_va,
        Y_aa,
        Y_fa,
        Y_c4,
        Y_pro,
        C_ch4,
        Y_ac,
        Y_h2,
        K_S_IN,
        k_m_su,
        K_S_su,
        pH_UL_aa,
        pH_LL_aa,
        k_m_aa,
        K_S_aa,
        k_m_fa,
        K_S_fa,
        K_I_h2_fa,
        k_m_c4,
        K_S_c4,
        K_I_h2_c4,
        k_m_pro,
        K_S_pro,
        K_I_h2_pro,
        k_m_ac,
        K_S_ac,
        K_I_nh3,
        pH_UL_ac,
        pH_LL_ac,
        k_m_h2,
        K_S_h2,
        pH_UL_h2,
        pH_LL_h2,
        k_dec_X_su,
        k_dec_X_aa,
        k_dec_X_fa,
        k_dec_X_c4,
        k_dec_X_pro,
        k_dec_X_ac,
        k_dec_X_h2,
        K_w,
        K_a_va,
        K_a_bu,
        K_a_pro,
        K_a_ac,
        K_a_co2,
        K_a_IN,
        k_A_B_va,
        k_A_B_bu,
        k_A_B_pro,
        k_A_B_ac,
        k_A_B_co2,
        k_A_B_IN,
        p_gas_h2o,
        k_p,
        k_L_a,
        K_H_co2,
        K_H_ch4,
        K_H_h2,
    )

pyadm1.core.adm_equations.BiochemicalProcesses

Combined biochemical process calculations for ADM1.

This class orchestrates the calculation of all process rates including inhibition factors and stoichiometric relationships.

Source code in pyadm1/core/adm_equations.py

class BiochemicalProcesses:
    """
    Combined biochemical process calculations for ADM1.

    This class orchestrates the calculation of all process rates including
    inhibition factors and stoichiometric relationships.
    """

    @staticmethod
    def calculate_inhibition_factors(
        S_H_ion: float,
        S_h2: float,
        S_nh4_ion: float,
        S_nh3: float,
        K_pH_aa: float,
        nn_aa: float,
        K_pH_ac: float,
        n_ac: float,
        K_pH_h2: float,
        n_h2: float,
        K_S_IN: float,
        K_I_h2_fa: float,
        K_I_h2_c4: float,
        K_I_h2_pro: float,
        K_I_nh3: float,
    ) -> Tuple[float, float, float, float, float, float, float, float, float]:
        """
        Calculate all inhibition factors for ADM1 processes.

        Args:
            S_H_ion: Hydrogen ion concentration [M]
            S_h2: Hydrogen gas concentration [kg COD/m³]
            S_nh4_ion: Ammonium concentration [M]
            S_nh3: Free ammonia concentration [M]
            K_pH_aa: pH inhibition constant for amino acid degraders [M]
            nn_aa: Hill coefficient for aa pH inhibition [-]
            K_pH_ac: pH inhibition constant for acetate degraders [M]
            n_ac: Hill coefficient for ac pH inhibition [-]
            K_pH_h2: pH inhibition constant for hydrogen degraders [M]
            n_h2: Hill coefficient for h2 pH inhibition [-]
            K_S_IN: Nitrogen half-saturation constant [M]
            K_I_h2_fa: H2 inhibition constant for LCFA degraders [kg COD/m³]
            K_I_h2_c4: H2 inhibition constant for C4 degraders [kg COD/m³]
            K_I_h2_pro: H2 inhibition constant for propionate degraders [kg COD/m³]
            K_I_nh3: Ammonia inhibition constant [M]

        Returns:
            Tuple of inhibition factors (I_pH_aa, I_pH_ac, I_pH_h2, I_IN_lim,
            I_h2_fa, I_h2_c4, I_h2_pro, I_nh3, I_5 through I_12)
        """
        inh = InhibitionFunctions()

        # pH inhibition factors
        I_pH_aa = inh.pH_inhibition(S_H_ion, K_pH_aa, nn_aa)
        I_pH_ac = inh.pH_inhibition(S_H_ion, K_pH_ac, n_ac)
        I_pH_h2 = inh.pH_inhibition(S_H_ion, K_pH_h2, n_h2)

        # Nitrogen limitation
        I_IN_lim = inh.nitrogen_limitation(S_nh4_ion, S_nh3, K_S_IN)

        # Hydrogen inhibition
        I_h2_fa = inh.hydrogen_inhibition(S_h2, K_I_h2_fa)
        I_h2_c4 = inh.hydrogen_inhibition(S_h2, K_I_h2_c4)
        I_h2_pro = inh.hydrogen_inhibition(S_h2, K_I_h2_pro)

        # Ammonia inhibition
        I_nh3 = inh.ammonia_inhibition(S_nh3, K_I_nh3)

        # Combined inhibition factors for different processes
        I_5 = I_pH_aa * I_IN_lim  # Sugar uptake
        I_6 = I_5  # Amino acid uptake
        I_7 = I_pH_aa * I_IN_lim * I_h2_fa  # LCFA uptake
        I_8 = I_pH_aa * I_IN_lim * I_h2_c4  # Valerate uptake
        I_9 = I_8  # Butyrate uptake
        I_10 = I_pH_aa * I_IN_lim * I_h2_pro  # Propionate uptake
        I_11 = I_pH_ac * I_IN_lim * I_nh3  # Acetate uptake
        I_12 = I_pH_h2 * I_IN_lim  # Hydrogen uptake

        return (
            I_pH_aa,
            I_pH_ac,
            I_pH_h2,
            I_IN_lim,
            I_h2_fa,
            I_h2_c4,
            I_h2_pro,
            I_nh3,
            I_5,
            I_6,
            I_7,
            I_8,
            I_9,
            I_10,
            I_11,
            I_12,
        )

    @staticmethod
    def calculate_process_rates(
        state: List[float],
        inhibitions: Tuple[float, ...],
        kinetic_params: dict,
        substrate_params: dict,
        hydro_factor: float = 1.0,
    ) -> Tuple[float, ...]:
        """
        Calculate all 19 biochemical process rates for ADM1.

        Args:
            state: ADM1 state vector (37 elements)
            inhibitions: Tuple of inhibition factors from calculate_inhibition_factors
            kinetic_params: Dictionary of kinetic parameters (k_m, K_S, k_dec, etc.)
            substrate_params: Dictionary of substrate-dependent parameters
                (k_dis, k_hyd_ch, k_hyd_pr, k_hyd_li)
            hydro_factor: Optional TS-dependent hydrolysis factor [-]. Values
                above 1 are interpreted as digester TS [%].

        Returns:
            Tuple of 19 process rates (Rho_1 through Rho_19)
        """
        # Unpack state
        S_su, S_aa, S_fa, S_va, S_bu, S_pro, S_ac, S_h2 = state[0:8]
        X_xc, X_ch, X_pr, X_li = state[12:16]
        X_su, X_aa, X_fa, X_c4, X_pro, X_ac, X_h2 = state[16:23]

        # Unpack inhibition factors (using indices 8-15 for combined factors)
        I_5, I_6, I_7, I_8, I_9, I_10, I_11, I_12 = inhibitions[8:16]

        # Unpack parameters
        k_dis = substrate_params["k_dis"]
        k_hyd_ch = substrate_params["k_hyd_ch"]
        k_hyd_pr = substrate_params["k_hyd_pr"]
        k_hyd_li = substrate_params["k_hyd_li"]

        if hydro_factor > 1.0:
            hydro_factor = 1.0 / (1.0 + (hydro_factor / 5.5) ** 2.3)
        elif hydro_factor < 0.0:
            hydro_factor = 0.0

        proc = ProcessRates()

        # Process rates (Rosen et al. 2006, BSM2)
        Rho_1 = proc.disintegration_rate(k_dis, X_xc)
        Rho_2 = proc.hydrolysis_rate(k_hyd_ch, X_ch, hydro_factor)
        Rho_3 = proc.hydrolysis_rate(k_hyd_pr, X_pr, hydro_factor)
        Rho_4 = proc.hydrolysis_rate(k_hyd_li, X_li, hydro_factor)

        # Uptake rates
        Rho_5 = proc.uptake_rate(kinetic_params["k_m_su"], S_su, kinetic_params["K_S_su"], X_su, I_5)
        Rho_6 = proc.uptake_rate(kinetic_params["k_m_aa"], S_aa, kinetic_params["K_S_aa"], X_aa, I_6)
        Rho_7 = proc.uptake_rate(kinetic_params["k_m_fa"], S_fa, kinetic_params["K_S_fa"], X_fa, I_7)

        # Valerate and butyrate uptake (with competition)
        competition_va = S_va / (S_bu + S_va + 1e-6)
        competition_bu = S_bu / (S_bu + S_va + 1e-6)
        Rho_8 = kinetic_params["k_m_c4"] * (S_va / (kinetic_params["K_S_c4"] + S_va)) * X_c4 * competition_va * I_8
        Rho_9 = kinetic_params["k_m_c4"] * (S_bu / (kinetic_params["K_S_c4"] + S_bu)) * X_c4 * competition_bu * I_9

        Rho_10 = proc.uptake_rate(kinetic_params["k_m_pro"], S_pro, kinetic_params["K_S_pro"], X_pro, I_10)
        Rho_11 = proc.uptake_rate(kinetic_params["k_m_ac"], S_ac, kinetic_params["K_S_ac"], X_ac, I_11)
        Rho_12 = proc.uptake_rate(kinetic_params["k_m_h2"], S_h2, kinetic_params["K_S_h2"], X_h2, I_12)

        # Decay rates
        Rho_13 = proc.decay_rate(kinetic_params["k_dec_X_su"], X_su)
        Rho_14 = proc.decay_rate(kinetic_params["k_dec_X_aa"], X_aa)
        Rho_15 = proc.decay_rate(kinetic_params["k_dec_X_fa"], X_fa)
        Rho_16 = proc.decay_rate(kinetic_params["k_dec_X_c4"], X_c4)
        Rho_17 = proc.decay_rate(kinetic_params["k_dec_X_pro"], X_pro)
        Rho_18 = proc.decay_rate(kinetic_params["k_dec_X_ac"], X_ac)
        Rho_19 = proc.decay_rate(kinetic_params["k_dec_X_h2"], X_h2)

        return (
            Rho_1,
            Rho_2,
            Rho_3,
            Rho_4,
            Rho_5,
            Rho_6,
            Rho_7,
            Rho_8,
            Rho_9,
            Rho_10,
            Rho_11,
            Rho_12,
            Rho_13,
            Rho_14,
            Rho_15,
            Rho_16,
            Rho_17,
            Rho_18,
            Rho_19,
        )

    @staticmethod
    def calculate_acid_base_rates(
        state: List[float], acid_base_params: dict
    ) -> Tuple[float, float, float, float, float, float]:
        """
        Calculate acid-base reaction rates for ODE implementation.

        Args:
            state: ADM1 state vector (37 elements)
            acid_base_params: Dictionary containing K_a and k_AB values

        Returns:
            Tuple of 6 acid-base rates (Rho_A_4 through Rho_A_11)
        """
        # Unpack state
        S_va, S_bu, S_pro, S_ac, S_co2, S_nh4_ion = (
            state[3],
            state[4],
            state[5],
            state[6],
            state[9],
            state[10],
        )
        S_va_ion, S_bu_ion, S_pro_ion, S_ac_ion = state[27:31]
        S_hco3_ion, S_nh3 = state[31], state[32]

        # Get pH to calculate H+ concentration
        # Note: In actual implementation, pH should be calculated from state
        # For now, we'll use a placeholder that should be replaced
        from biogas import ADMstate

        S_H_ion = 10 ** (-ADMstate.calcPHOfADMstate(state))

        ab = AcidBaseKinetics()

        # Acid-base rates
        Rho_A_4 = ab.acid_base_rate(
            acid_base_params["k_A_B_va"],
            S_va_ion,
            S_H_ion,
            acid_base_params["K_a_va"],
            S_va - S_va_ion,
        )

        Rho_A_5 = ab.acid_base_rate(
            acid_base_params["k_A_B_bu"],
            S_bu_ion,
            S_H_ion,
            acid_base_params["K_a_bu"],
            S_bu - S_bu_ion,
        )

        Rho_A_6 = ab.acid_base_rate(
            acid_base_params["k_A_B_pro"],
            S_pro_ion,
            S_H_ion,
            acid_base_params["K_a_pro"],
            S_pro - S_pro_ion,
        )

        Rho_A_7 = ab.acid_base_rate(
            acid_base_params["k_A_B_ac"],
            S_ac_ion,
            S_H_ion,
            acid_base_params["K_a_ac"],
            S_ac - S_ac_ion,
        )

        Rho_A_10 = ab.acid_base_rate(
            acid_base_params["k_A_B_co2"],
            S_hco3_ion,
            S_H_ion,
            acid_base_params["K_a_co2"],
            S_co2,
        )

        Rho_A_11 = ab.acid_base_rate(
            acid_base_params["k_A_B_IN"],
            S_nh3,
            S_H_ion,
            acid_base_params["K_a_IN"],
            S_nh4_ion,
        )

        return (Rho_A_4, Rho_A_5, Rho_A_6, Rho_A_7, Rho_A_10, Rho_A_11)

    @staticmethod
    def calculate_gas_transfer_rates(
        state: List[float],
        gas_params: dict,
        RT: float,
        V_liq: float,
        V_gas: float,
        p_ext: float,
    ) -> Tuple[float, float, float, float]:
        """
        Calculate gas-liquid transfer and gas outlet rates.

        Args:
            state: ADM1 state vector (37 elements)
            gas_params: Dictionary containing k_L_a, K_H constants, k_p
            RT: Gas constant × temperature [bar·m³/kmol]
            V_liq: Liquid volume [m³]
            V_gas: Gas volume [m³]
            p_ext: External pressure [bar]
        Returns:
            Tuple of 4 rates (Rho_T_8, Rho_T_9, Rho_T_10, Rho_T_11)
        """
        # Unpack state
        S_h2, S_ch4, S_co2 = state[7:10]
        p_gas_h2, p_gas_ch4, p_gas_co2, pTOTAL = state[33:37]

        gt = GasTransfer()

        # Gas transfer rates
        Rho_T_8 = gt.gas_transfer_rate(
            gas_params["k_L_a"],
            S_h2,
            p_gas_h2,
            gas_params["K_H_h2"],
            RT,
            16.0,
            V_liq,
            V_gas,
        )

        Rho_T_9 = gt.gas_transfer_rate(
            gas_params["k_L_a"],
            S_ch4,
            p_gas_ch4,
            gas_params["K_H_ch4"],
            RT,
            64.0,
            V_liq,
            V_gas,
        )

        Rho_T_10 = gt.gas_transfer_rate(
            gas_params["k_L_a"],
            S_co2,
            p_gas_co2,
            gas_params["K_H_co2"],
            RT,
            1.0,
            V_liq,
            V_gas,
        )

        # Gas outlet rate
        Rho_T_11 = gt.gas_outlet_rate(gas_params["k_p"], pTOTAL, p_ext, V_liq, V_gas)

        return (Rho_T_8, Rho_T_9, Rho_T_10, Rho_T_11)

Functions

calculate_acid_base_rates(state, acid_base_params) staticmethod

Calculate acid-base reaction rates for ODE implementation.

Parameters:

Name	Type	Description	Default
state	List[float]	ADM1 state vector (37 elements)	required
acid_base_params	dict	Dictionary containing K_a and k_AB values	required

Returns:

Type	Description
Tuple[float, float, float, float, float, float]	Tuple of 6 acid-base rates (Rho_A_4 through Rho_A_11)

Source code in pyadm1/core/adm_equations.py

@staticmethod
def calculate_acid_base_rates(
    state: List[float], acid_base_params: dict
) -> Tuple[float, float, float, float, float, float]:
    """
    Calculate acid-base reaction rates for ODE implementation.

    Args:
        state: ADM1 state vector (37 elements)
        acid_base_params: Dictionary containing K_a and k_AB values

    Returns:
        Tuple of 6 acid-base rates (Rho_A_4 through Rho_A_11)
    """
    # Unpack state
    S_va, S_bu, S_pro, S_ac, S_co2, S_nh4_ion = (
        state[3],
        state[4],
        state[5],
        state[6],
        state[9],
        state[10],
    )
    S_va_ion, S_bu_ion, S_pro_ion, S_ac_ion = state[27:31]
    S_hco3_ion, S_nh3 = state[31], state[32]

    # Get pH to calculate H+ concentration
    # Note: In actual implementation, pH should be calculated from state
    # For now, we'll use a placeholder that should be replaced
    from biogas import ADMstate

    S_H_ion = 10 ** (-ADMstate.calcPHOfADMstate(state))

    ab = AcidBaseKinetics()

    # Acid-base rates
    Rho_A_4 = ab.acid_base_rate(
        acid_base_params["k_A_B_va"],
        S_va_ion,
        S_H_ion,
        acid_base_params["K_a_va"],
        S_va - S_va_ion,
    )

    Rho_A_5 = ab.acid_base_rate(
        acid_base_params["k_A_B_bu"],
        S_bu_ion,
        S_H_ion,
        acid_base_params["K_a_bu"],
        S_bu - S_bu_ion,
    )

    Rho_A_6 = ab.acid_base_rate(
        acid_base_params["k_A_B_pro"],
        S_pro_ion,
        S_H_ion,
        acid_base_params["K_a_pro"],
        S_pro - S_pro_ion,
    )

    Rho_A_7 = ab.acid_base_rate(
        acid_base_params["k_A_B_ac"],
        S_ac_ion,
        S_H_ion,
        acid_base_params["K_a_ac"],
        S_ac - S_ac_ion,
    )

    Rho_A_10 = ab.acid_base_rate(
        acid_base_params["k_A_B_co2"],
        S_hco3_ion,
        S_H_ion,
        acid_base_params["K_a_co2"],
        S_co2,
    )

    Rho_A_11 = ab.acid_base_rate(
        acid_base_params["k_A_B_IN"],
        S_nh3,
        S_H_ion,
        acid_base_params["K_a_IN"],
        S_nh4_ion,
    )

    return (Rho_A_4, Rho_A_5, Rho_A_6, Rho_A_7, Rho_A_10, Rho_A_11)

calculate_gas_transfer_rates(state, gas_params, RT, V_liq, V_gas, p_ext) staticmethod

Calculate gas-liquid transfer and gas outlet rates.

Parameters:

Name	Type	Description	Default
state	List[float]	ADM1 state vector (37 elements)	required
gas_params	dict	Dictionary containing k_L_a, K_H constants, k_p	required
RT	float	Gas constant × temperature [bar·m³/kmol]	required
V_liq	float	Liquid volume [m³]	required
V_gas	float	Gas volume [m³]	required
p_ext	float	External pressure [bar]	required

Returns: Tuple of 4 rates (Rho_T_8, Rho_T_9, Rho_T_10, Rho_T_11)

Source code in pyadm1/core/adm_equations.py

@staticmethod
def calculate_gas_transfer_rates(
    state: List[float],
    gas_params: dict,
    RT: float,
    V_liq: float,
    V_gas: float,
    p_ext: float,
) -> Tuple[float, float, float, float]:
    """
    Calculate gas-liquid transfer and gas outlet rates.

    Args:
        state: ADM1 state vector (37 elements)
        gas_params: Dictionary containing k_L_a, K_H constants, k_p
        RT: Gas constant × temperature [bar·m³/kmol]
        V_liq: Liquid volume [m³]
        V_gas: Gas volume [m³]
        p_ext: External pressure [bar]
    Returns:
        Tuple of 4 rates (Rho_T_8, Rho_T_9, Rho_T_10, Rho_T_11)
    """
    # Unpack state
    S_h2, S_ch4, S_co2 = state[7:10]
    p_gas_h2, p_gas_ch4, p_gas_co2, pTOTAL = state[33:37]

    gt = GasTransfer()

    # Gas transfer rates
    Rho_T_8 = gt.gas_transfer_rate(
        gas_params["k_L_a"],
        S_h2,
        p_gas_h2,
        gas_params["K_H_h2"],
        RT,
        16.0,
        V_liq,
        V_gas,
    )

    Rho_T_9 = gt.gas_transfer_rate(
        gas_params["k_L_a"],
        S_ch4,
        p_gas_ch4,
        gas_params["K_H_ch4"],
        RT,
        64.0,
        V_liq,
        V_gas,
    )

    Rho_T_10 = gt.gas_transfer_rate(
        gas_params["k_L_a"],
        S_co2,
        p_gas_co2,
        gas_params["K_H_co2"],
        RT,
        1.0,
        V_liq,
        V_gas,
    )

    # Gas outlet rate
    Rho_T_11 = gt.gas_outlet_rate(gas_params["k_p"], pTOTAL, p_ext, V_liq, V_gas)

    return (Rho_T_8, Rho_T_9, Rho_T_10, Rho_T_11)

calculate_inhibition_factors(S_H_ion, S_h2, S_nh4_ion, S_nh3, K_pH_aa, nn_aa, K_pH_ac, n_ac, K_pH_h2, n_h2, K_S_IN, K_I_h2_fa, K_I_h2_c4, K_I_h2_pro, K_I_nh3) staticmethod

Calculate all inhibition factors for ADM1 processes.

Parameters:

Name	Type	Description	Default
S_H_ion	float	Hydrogen ion concentration [M]	required
S_h2	float	Hydrogen gas concentration [kg COD/m³]	required
S_nh4_ion	float	Ammonium concentration [M]	required
S_nh3	float	Free ammonia concentration [M]	required
K_pH_aa	float	pH inhibition constant for amino acid degraders [M]	required
nn_aa	float	Hill coefficient for aa pH inhibition [-]	required
K_pH_ac	float	pH inhibition constant for acetate degraders [M]	required
n_ac	float	Hill coefficient for ac pH inhibition [-]	required
K_pH_h2	float	pH inhibition constant for hydrogen degraders [M]	required
n_h2	float	Hill coefficient for h2 pH inhibition [-]	required
K_S_IN	float	Nitrogen half-saturation constant [M]	required
K_I_h2_fa	float	H2 inhibition constant for LCFA degraders [kg COD/m³]	required
K_I_h2_c4	float	H2 inhibition constant for C4 degraders [kg COD/m³]	required
K_I_h2_pro	float	H2 inhibition constant for propionate degraders [kg COD/m³]	required
K_I_nh3	float	Ammonia inhibition constant [M]	required

Returns:

Type	Description
float	Tuple of inhibition factors (I_pH_aa, I_pH_ac, I_pH_h2, I_IN_lim,
float	I_h2_fa, I_h2_c4, I_h2_pro, I_nh3, I_5 through I_12)

Source code in pyadm1/core/adm_equations.py

@staticmethod
def calculate_inhibition_factors(
    S_H_ion: float,
    S_h2: float,
    S_nh4_ion: float,
    S_nh3: float,
    K_pH_aa: float,
    nn_aa: float,
    K_pH_ac: float,
    n_ac: float,
    K_pH_h2: float,
    n_h2: float,
    K_S_IN: float,
    K_I_h2_fa: float,
    K_I_h2_c4: float,
    K_I_h2_pro: float,
    K_I_nh3: float,
) -> Tuple[float, float, float, float, float, float, float, float, float]:
    """
    Calculate all inhibition factors for ADM1 processes.

    Args:
        S_H_ion: Hydrogen ion concentration [M]
        S_h2: Hydrogen gas concentration [kg COD/m³]
        S_nh4_ion: Ammonium concentration [M]
        S_nh3: Free ammonia concentration [M]
        K_pH_aa: pH inhibition constant for amino acid degraders [M]
        nn_aa: Hill coefficient for aa pH inhibition [-]
        K_pH_ac: pH inhibition constant for acetate degraders [M]
        n_ac: Hill coefficient for ac pH inhibition [-]
        K_pH_h2: pH inhibition constant for hydrogen degraders [M]
        n_h2: Hill coefficient for h2 pH inhibition [-]
        K_S_IN: Nitrogen half-saturation constant [M]
        K_I_h2_fa: H2 inhibition constant for LCFA degraders [kg COD/m³]
        K_I_h2_c4: H2 inhibition constant for C4 degraders [kg COD/m³]
        K_I_h2_pro: H2 inhibition constant for propionate degraders [kg COD/m³]
        K_I_nh3: Ammonia inhibition constant [M]

    Returns:
        Tuple of inhibition factors (I_pH_aa, I_pH_ac, I_pH_h2, I_IN_lim,
        I_h2_fa, I_h2_c4, I_h2_pro, I_nh3, I_5 through I_12)
    """
    inh = InhibitionFunctions()

    # pH inhibition factors
    I_pH_aa = inh.pH_inhibition(S_H_ion, K_pH_aa, nn_aa)
    I_pH_ac = inh.pH_inhibition(S_H_ion, K_pH_ac, n_ac)
    I_pH_h2 = inh.pH_inhibition(S_H_ion, K_pH_h2, n_h2)

    # Nitrogen limitation
    I_IN_lim = inh.nitrogen_limitation(S_nh4_ion, S_nh3, K_S_IN)

    # Hydrogen inhibition
    I_h2_fa = inh.hydrogen_inhibition(S_h2, K_I_h2_fa)
    I_h2_c4 = inh.hydrogen_inhibition(S_h2, K_I_h2_c4)
    I_h2_pro = inh.hydrogen_inhibition(S_h2, K_I_h2_pro)

    # Ammonia inhibition
    I_nh3 = inh.ammonia_inhibition(S_nh3, K_I_nh3)

    # Combined inhibition factors for different processes
    I_5 = I_pH_aa * I_IN_lim  # Sugar uptake
    I_6 = I_5  # Amino acid uptake
    I_7 = I_pH_aa * I_IN_lim * I_h2_fa  # LCFA uptake
    I_8 = I_pH_aa * I_IN_lim * I_h2_c4  # Valerate uptake
    I_9 = I_8  # Butyrate uptake
    I_10 = I_pH_aa * I_IN_lim * I_h2_pro  # Propionate uptake
    I_11 = I_pH_ac * I_IN_lim * I_nh3  # Acetate uptake
    I_12 = I_pH_h2 * I_IN_lim  # Hydrogen uptake

    return (
        I_pH_aa,
        I_pH_ac,
        I_pH_h2,
        I_IN_lim,
        I_h2_fa,
        I_h2_c4,
        I_h2_pro,
        I_nh3,
        I_5,
        I_6,
        I_7,
        I_8,
        I_9,
        I_10,
        I_11,
        I_12,
    )

calculate_process_rates(state, inhibitions, kinetic_params, substrate_params, hydro_factor=1.0) staticmethod

Calculate all 19 biochemical process rates for ADM1.

Parameters:

Name	Type	Description	Default
state	List[float]	ADM1 state vector (37 elements)	required
inhibitions	Tuple[float, ...]	Tuple of inhibition factors from calculate_inhibition_factors	required
kinetic_params	dict	Dictionary of kinetic parameters (k_m, K_S, k_dec, etc.)	required
substrate_params	dict	Dictionary of substrate-dependent parameters (k_dis, k_hyd_ch, k_hyd_pr, k_hyd_li)	required
hydro_factor	float	Optional TS-dependent hydrolysis factor [-]. Values above 1 are interpreted as digester TS [%].	1.0

Returns:

Type	Description
Tuple[float, ...]	Tuple of 19 process rates (Rho_1 through Rho_19)

Source code in pyadm1/core/adm_equations.py

@staticmethod
def calculate_process_rates(
    state: List[float],
    inhibitions: Tuple[float, ...],
    kinetic_params: dict,
    substrate_params: dict,
    hydro_factor: float = 1.0,
) -> Tuple[float, ...]:
    """
    Calculate all 19 biochemical process rates for ADM1.

    Args:
        state: ADM1 state vector (37 elements)
        inhibitions: Tuple of inhibition factors from calculate_inhibition_factors
        kinetic_params: Dictionary of kinetic parameters (k_m, K_S, k_dec, etc.)
        substrate_params: Dictionary of substrate-dependent parameters
            (k_dis, k_hyd_ch, k_hyd_pr, k_hyd_li)
        hydro_factor: Optional TS-dependent hydrolysis factor [-]. Values
            above 1 are interpreted as digester TS [%].

    Returns:
        Tuple of 19 process rates (Rho_1 through Rho_19)
    """
    # Unpack state
    S_su, S_aa, S_fa, S_va, S_bu, S_pro, S_ac, S_h2 = state[0:8]
    X_xc, X_ch, X_pr, X_li = state[12:16]
    X_su, X_aa, X_fa, X_c4, X_pro, X_ac, X_h2 = state[16:23]

    # Unpack inhibition factors (using indices 8-15 for combined factors)
    I_5, I_6, I_7, I_8, I_9, I_10, I_11, I_12 = inhibitions[8:16]

    # Unpack parameters
    k_dis = substrate_params["k_dis"]
    k_hyd_ch = substrate_params["k_hyd_ch"]
    k_hyd_pr = substrate_params["k_hyd_pr"]
    k_hyd_li = substrate_params["k_hyd_li"]

    if hydro_factor > 1.0:
        hydro_factor = 1.0 / (1.0 + (hydro_factor / 5.5) ** 2.3)
    elif hydro_factor < 0.0:
        hydro_factor = 0.0

    proc = ProcessRates()

    # Process rates (Rosen et al. 2006, BSM2)
    Rho_1 = proc.disintegration_rate(k_dis, X_xc)
    Rho_2 = proc.hydrolysis_rate(k_hyd_ch, X_ch, hydro_factor)
    Rho_3 = proc.hydrolysis_rate(k_hyd_pr, X_pr, hydro_factor)
    Rho_4 = proc.hydrolysis_rate(k_hyd_li, X_li, hydro_factor)

    # Uptake rates
    Rho_5 = proc.uptake_rate(kinetic_params["k_m_su"], S_su, kinetic_params["K_S_su"], X_su, I_5)
    Rho_6 = proc.uptake_rate(kinetic_params["k_m_aa"], S_aa, kinetic_params["K_S_aa"], X_aa, I_6)
    Rho_7 = proc.uptake_rate(kinetic_params["k_m_fa"], S_fa, kinetic_params["K_S_fa"], X_fa, I_7)

    # Valerate and butyrate uptake (with competition)
    competition_va = S_va / (S_bu + S_va + 1e-6)
    competition_bu = S_bu / (S_bu + S_va + 1e-6)
    Rho_8 = kinetic_params["k_m_c4"] * (S_va / (kinetic_params["K_S_c4"] + S_va)) * X_c4 * competition_va * I_8
    Rho_9 = kinetic_params["k_m_c4"] * (S_bu / (kinetic_params["K_S_c4"] + S_bu)) * X_c4 * competition_bu * I_9

    Rho_10 = proc.uptake_rate(kinetic_params["k_m_pro"], S_pro, kinetic_params["K_S_pro"], X_pro, I_10)
    Rho_11 = proc.uptake_rate(kinetic_params["k_m_ac"], S_ac, kinetic_params["K_S_ac"], X_ac, I_11)
    Rho_12 = proc.uptake_rate(kinetic_params["k_m_h2"], S_h2, kinetic_params["K_S_h2"], X_h2, I_12)

    # Decay rates
    Rho_13 = proc.decay_rate(kinetic_params["k_dec_X_su"], X_su)
    Rho_14 = proc.decay_rate(kinetic_params["k_dec_X_aa"], X_aa)
    Rho_15 = proc.decay_rate(kinetic_params["k_dec_X_fa"], X_fa)
    Rho_16 = proc.decay_rate(kinetic_params["k_dec_X_c4"], X_c4)
    Rho_17 = proc.decay_rate(kinetic_params["k_dec_X_pro"], X_pro)
    Rho_18 = proc.decay_rate(kinetic_params["k_dec_X_ac"], X_ac)
    Rho_19 = proc.decay_rate(kinetic_params["k_dec_X_h2"], X_h2)

    return (
        Rho_1,
        Rho_2,
        Rho_3,
        Rho_4,
        Rho_5,
        Rho_6,
        Rho_7,
        Rho_8,
        Rho_9,
        Rho_10,
        Rho_11,
        Rho_12,
        Rho_13,
        Rho_14,
        Rho_15,
        Rho_16,
        Rho_17,
        Rho_18,
        Rho_19,
    )

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