Simulation API

pyadm1.simulation.simulator.Simulator

Single-plant ADM1 simulation driver.

Attributes

adm1 : ADM1 The model instance being driven. solver : ODESolver ODE solver wrapper (default BDF, rtol=1e-6, atol=1e-8).

Source code in pyadm1/simulation/simulator.py

class Simulator:
    """
    Single-plant ADM1 simulation driver.

    Attributes
    ----------
    adm1 : ADM1
        The model instance being driven.
    solver : ODESolver
        ODE solver wrapper (default BDF, rtol=1e-6, atol=1e-8).
    """

    def __init__(self, adm1: ADM1, solver: ODESolver = None) -> None:
        self._adm1 = adm1
        self._solver = solver or create_solver(method="BDF", rtol=1e-6, atol=1e-8)

    def simulate_AD_plant(self, tstep: List[float], state_zero: List[float]) -> List[float]:
        """
        Integrate the ADM1 ODE for the requested time span.

        Parameters
        ----------
        tstep : [t_start, t_end]
            Time span [days].
        state_zero : list of float
            Initial state vector (41 elements).

        Returns
        -------
        list of float
            Final state vector after integration.
        """
        result = self._solver.solve(fun=self._adm1.ADM_ODE, t_span=tstep, y0=state_zero)
        final_state = result.y[:, -1].tolist()

        # Update tracked process indicators (pH, gas) so that downstream code
        # can read them from the ADM1 instance directly.
        self._adm1.print_params_at_current_state(final_state)

        return final_state

    def simulate_gas_production(
        self,
        tstep: List[float],
        state_zero: List[float],
        Q: List[float],
    ) -> Tuple[float, float]:
        """
        Integrate the ODE and return final-state gas production rates.

        Parameters
        ----------
        tstep : [t_start, t_end]
            Time span [days].
        state_zero : list of float
            Initial state vector (41 elements).
        Q : list of float
            Substrate feed rates [m³/d].

        Returns
        -------
        (q_gas, q_ch4) : tuple of float
            Total biogas flow and methane flow [Nm³/d].
        """
        self._adm1.create_influent(Q, 0)
        result = self._solver.solve(fun=self._adm1.ADM_ODE, t_span=tstep, y0=state_zero)
        final_state = result.y[:, -1]
        q_gas, q_ch4, _, _, _ = self._adm1.calc_gas(
            float(final_state[37]),
            float(final_state[38]),
            float(final_state[39]),
            float(final_state[40]),
        )
        return float(q_gas), float(q_ch4)

    @property
    def adm1(self) -> ADM1:
        """ADM1 model instance."""
        return self._adm1

    @property
    def solver(self) -> ODESolver:
        """ODE solver instance."""
        return self._solver

Attributes

adm1 property

ADM1 model instance.

solver property

ODE solver instance.

Functions

simulate_AD_plant(tstep, state_zero)

Integrate the ADM1 ODE for the requested time span.

Parameters

tstep : [t_start, t_end] Time span [days]. state_zero : list of float Initial state vector (41 elements).

Returns

list of float Final state vector after integration.

Source code in pyadm1/simulation/simulator.py

def simulate_AD_plant(self, tstep: List[float], state_zero: List[float]) -> List[float]:
    """
    Integrate the ADM1 ODE for the requested time span.

    Parameters
    ----------
    tstep : [t_start, t_end]
        Time span [days].
    state_zero : list of float
        Initial state vector (41 elements).

    Returns
    -------
    list of float
        Final state vector after integration.
    """
    result = self._solver.solve(fun=self._adm1.ADM_ODE, t_span=tstep, y0=state_zero)
    final_state = result.y[:, -1].tolist()

    # Update tracked process indicators (pH, gas) so that downstream code
    # can read them from the ADM1 instance directly.
    self._adm1.print_params_at_current_state(final_state)

    return final_state

simulate_gas_production(tstep, state_zero, Q)

Integrate the ODE and return final-state gas production rates.

Parameters

tstep : [t_start, t_end] Time span [days]. state_zero : list of float Initial state vector (41 elements). Q : list of float Substrate feed rates [m³/d].

Returns

(q_gas, q_ch4) : tuple of float Total biogas flow and methane flow [Nm³/d].

Source code in pyadm1/simulation/simulator.py

def simulate_gas_production(
    self,
    tstep: List[float],
    state_zero: List[float],
    Q: List[float],
) -> Tuple[float, float]:
    """
    Integrate the ODE and return final-state gas production rates.

    Parameters
    ----------
    tstep : [t_start, t_end]
        Time span [days].
    state_zero : list of float
        Initial state vector (41 elements).
    Q : list of float
        Substrate feed rates [m³/d].

    Returns
    -------
    (q_gas, q_ch4) : tuple of float
        Total biogas flow and methane flow [Nm³/d].
    """
    self._adm1.create_influent(Q, 0)
    result = self._solver.solve(fun=self._adm1.ADM_ODE, t_span=tstep, y0=state_zero)
    final_state = result.y[:, -1]
    q_gas, q_ch4, _, _, _ = self._adm1.calc_gas(
        float(final_state[37]),
        float(final_state[38]),
        float(final_state[39]),
        float(final_state[40]),
    )
    return float(q_gas), float(q_ch4)

pyadm1.simulation.parallel.ParallelSimulator

Parallel simulator for running multiple ADM1 scenarios concurrently.

Source code in pyadm1/simulation/parallel.py

class ParallelSimulator:
    """
    Parallel simulator for running multiple ADM1 scenarios concurrently.
    """

    def __init__(self, adm1: "ADM1", n_workers: Optional[int] = None, verbose: bool = True):
        """
        Parameters
        ----------
        adm1 : ADM1
            Base model instance — its V_liq, V_gas, T_ad and feedstock
            substrate IDs are copied to each worker.
        n_workers : int, optional
            Worker process count (default = ``cpu_count() - 1``).
        verbose : bool
            Whether to print progress.
        """
        self.adm1 = adm1
        self.n_workers = n_workers or max(1, mp.cpu_count() - 1)
        self.verbose = verbose

    def run_scenarios(
        self,
        scenarios: List[Dict[str, Any]],
        duration: float,
        initial_state: List[float],
        dt: float = 1.0 / 24.0,
        compute_metrics: bool = True,
        save_time_series: bool = False,
    ) -> List[ScenarioResult]:
        """Run multiple scenarios in parallel."""
        if self.verbose:
            print(f"Starting parallel simulation with {len(scenarios)} scenarios")
            print(f"Using {self.n_workers} worker processes")

        start_time = time.time()

        worker_func = partial(
            _run_single_scenario,
            adm1_config=self._serialize_adm1(),
            duration=duration,
            initial_state=initial_state,
            dt=dt,
            compute_metrics=compute_metrics,
            save_time_series=save_time_series,
            verbose=self.verbose,
        )

        scenarios_with_ids = [(i, scenario) for i, scenario in enumerate(scenarios)]

        use_sequential = self.n_workers <= 1 or len(scenarios_with_ids) <= 1

        if use_sequential:
            results = []
            for i, scenario_data in enumerate(scenarios_with_ids):
                results.append(worker_func(scenario_data))
                if self.verbose and ((i + 1) % 10 == 0 or (i + 1) == len(scenarios)):
                    print(f"  Completed {i + 1}/{len(scenarios)} scenarios")
        else:
            ctx = _get_mp_context()
            with ctx.Pool(processes=self.n_workers) as pool:
                if self.verbose:
                    results = []
                    for i, result in enumerate(pool.imap(worker_func, scenarios_with_ids)):
                        results.append(result)
                        if (i + 1) % 10 == 0 or (i + 1) == len(scenarios):
                            print(f"  Completed {i + 1}/{len(scenarios)} scenarios")
                else:
                    results = pool.map(worker_func, scenarios_with_ids)

        elapsed_time = time.time() - start_time

        if self.verbose:
            n_success = sum(1 for r in results if r.success)
            print("\nSimulation complete:")
            print(f"  Total scenarios: {len(scenarios)}")
            print(f"  Successful: {n_success}")
            print(f"  Failed: {len(scenarios) - n_success}")
            print(f"  Total time: {elapsed_time:.1f} seconds")
            print(f"  Average time per scenario: {elapsed_time / len(scenarios):.2f} seconds")

        return results

    def parameter_sweep(
        self,
        config: ParameterSweepConfig,
        duration: float,
        initial_state: List[float],
        **kwargs: Any,
    ) -> List[ScenarioResult]:
        """Run a single-parameter sweep."""
        scenarios = []
        for value in config.values:
            scenario = config.other_params.copy()
            scenario[config.parameter_name] = value
            scenarios.append(scenario)

        if self.verbose:
            print(f"Parameter sweep: {config.parameter_name}")
            print(f"  Values: {config.values}")

        return self.run_scenarios(scenarios, duration, initial_state, **kwargs)

    def multi_parameter_sweep(
        self,
        parameter_configs: Dict[str, List[float]],
        duration: float,
        initial_state: List[float],
        fixed_params: Optional[Dict[str, Any]] = None,
        **kwargs: Any,
    ) -> List[ScenarioResult]:
        """Run a multi-parameter sweep (full factorial design)."""
        fixed_params = fixed_params or {}

        param_names = list(parameter_configs.keys())
        param_values = [parameter_configs[name] for name in param_names]

        scenarios = []
        for combination in self._generate_combinations(param_values):
            scenario = fixed_params.copy()
            for name, value in zip(param_names, combination):
                scenario[name] = value
            scenarios.append(scenario)

        if self.verbose:
            print("Multi-parameter sweep:")
            for name, values in parameter_configs.items():
                print(f"  {name}: {len(values)} values")
            print(f"  Total combinations: {len(scenarios)}")

        return self.run_scenarios(scenarios, duration, initial_state, **kwargs)

    def monte_carlo(
        self,
        config: MonteCarloConfig,
        duration: float,
        initial_state: List[float],
        **kwargs: Any,
    ) -> List[ScenarioResult]:
        """Run Monte Carlo simulation with parameter uncertainty."""
        if config.seed is not None:
            np.random.seed(config.seed)

        scenarios = []
        for _ in range(config.n_samples):
            scenario = config.fixed_params.copy()
            for param_name, (mean, std) in config.parameter_distributions.items():
                value = np.random.normal(mean, std)
                scenario[param_name] = value
            scenarios.append(scenario)

        if self.verbose:
            print("Monte Carlo simulation:")
            print(f"  Samples: {config.n_samples}")
            print("  Parameters with uncertainty:")
            for param_name, (mean, std) in config.parameter_distributions.items():
                print(f"    {param_name}: N({mean}, {std}²)")

        return self.run_scenarios(scenarios, duration, initial_state, **kwargs)

    def summarize_results(self, results: List[ScenarioResult], metrics: Optional[List[str]] = None) -> Dict[str, Any]:
        """Compute summary statistics across multiple scenarios."""
        successful = [r for r in results if r.success]
        n_scenarios = len(results)
        n_successful = len(successful)

        summary = {
            "n_scenarios": n_scenarios,
            "n_successful": n_successful,
            "n_failed": n_scenarios - n_successful,
            "success_rate": n_successful / n_scenarios if n_scenarios > 0 else 0,
            "metrics": {},
        }

        if not successful:
            if n_scenarios > 0:
                summary["error"] = "No successful scenarios"
            else:
                summary["error"] = "No scenarios to summarize"
            return summary

        if metrics is None:
            metrics = list(successful[0].metrics.keys())

        for metric_name in metrics:
            values = [r.metrics.get(metric_name, np.nan) for r in successful]
            values = [v for v in values if not np.isnan(v)]

            if values:
                summary["metrics"][metric_name] = {
                    "mean": np.mean(values),
                    "std": np.std(values),
                    "min": np.min(values),
                    "max": np.max(values),
                    "median": np.median(values),
                    "q25": np.percentile(values, 25),
                    "q75": np.percentile(values, 75),
                }

        return summary

    def _serialize_adm1(self) -> Dict[str, Any]:
        """Serialize ADM1 model configuration for the worker pool."""
        feedstock = self.adm1.feedstock
        substrate_ids: List[str] = []
        feeding_freq = 24
        if feedstock is not None:
            substrate_ids = list(getattr(feedstock, "substrate_ids", []))
            feeding_freq = int(getattr(feedstock, "feeding_freq", 24))

        return {
            "V_liq": self.adm1.V_liq,
            "V_gas": self.adm1._V_gas,
            "T_ad": self.adm1._T_ad,
            "feedstock_substrates": substrate_ids,
            "feeding_freq": feeding_freq,
        }

    @staticmethod
    def _generate_combinations(value_lists: List[List[Any]]) -> List[Tuple]:
        """Cartesian product of a list of value lists."""
        if not value_lists:
            return [()]

        result = []
        for value in value_lists[0]:
            for rest in ParallelSimulator._generate_combinations(value_lists[1:]):
                result.append((value,) + rest)

        return result

Functions

__init__(adm1, n_workers=None, verbose=True)

Parameters

adm1 : ADM1 Base model instance — its V_liq, V_gas, T_ad and feedstock substrate IDs are copied to each worker. n_workers : int, optional Worker process count (default = cpu_count() - 1). verbose : bool Whether to print progress.

Source code in pyadm1/simulation/parallel.py

def __init__(self, adm1: "ADM1", n_workers: Optional[int] = None, verbose: bool = True):
    """
    Parameters
    ----------
    adm1 : ADM1
        Base model instance — its V_liq, V_gas, T_ad and feedstock
        substrate IDs are copied to each worker.
    n_workers : int, optional
        Worker process count (default = ``cpu_count() - 1``).
    verbose : bool
        Whether to print progress.
    """
    self.adm1 = adm1
    self.n_workers = n_workers or max(1, mp.cpu_count() - 1)
    self.verbose = verbose

monte_carlo(config, duration, initial_state, **kwargs)

Run Monte Carlo simulation with parameter uncertainty.

Source code in pyadm1/simulation/parallel.py

def monte_carlo(
    self,
    config: MonteCarloConfig,
    duration: float,
    initial_state: List[float],
    **kwargs: Any,
) -> List[ScenarioResult]:
    """Run Monte Carlo simulation with parameter uncertainty."""
    if config.seed is not None:
        np.random.seed(config.seed)

    scenarios = []
    for _ in range(config.n_samples):
        scenario = config.fixed_params.copy()
        for param_name, (mean, std) in config.parameter_distributions.items():
            value = np.random.normal(mean, std)
            scenario[param_name] = value
        scenarios.append(scenario)

    if self.verbose:
        print("Monte Carlo simulation:")
        print(f"  Samples: {config.n_samples}")
        print("  Parameters with uncertainty:")
        for param_name, (mean, std) in config.parameter_distributions.items():
            print(f"    {param_name}: N({mean}, {std}²)")

    return self.run_scenarios(scenarios, duration, initial_state, **kwargs)

multi_parameter_sweep(parameter_configs, duration, initial_state, fixed_params=None, **kwargs)

Run a multi-parameter sweep (full factorial design).

Source code in pyadm1/simulation/parallel.py

def multi_parameter_sweep(
    self,
    parameter_configs: Dict[str, List[float]],
    duration: float,
    initial_state: List[float],
    fixed_params: Optional[Dict[str, Any]] = None,
    **kwargs: Any,
) -> List[ScenarioResult]:
    """Run a multi-parameter sweep (full factorial design)."""
    fixed_params = fixed_params or {}

    param_names = list(parameter_configs.keys())
    param_values = [parameter_configs[name] for name in param_names]

    scenarios = []
    for combination in self._generate_combinations(param_values):
        scenario = fixed_params.copy()
        for name, value in zip(param_names, combination):
            scenario[name] = value
        scenarios.append(scenario)

    if self.verbose:
        print("Multi-parameter sweep:")
        for name, values in parameter_configs.items():
            print(f"  {name}: {len(values)} values")
        print(f"  Total combinations: {len(scenarios)}")

    return self.run_scenarios(scenarios, duration, initial_state, **kwargs)

parameter_sweep(config, duration, initial_state, **kwargs)

Run a single-parameter sweep.

Source code in pyadm1/simulation/parallel.py

def parameter_sweep(
    self,
    config: ParameterSweepConfig,
    duration: float,
    initial_state: List[float],
    **kwargs: Any,
) -> List[ScenarioResult]:
    """Run a single-parameter sweep."""
    scenarios = []
    for value in config.values:
        scenario = config.other_params.copy()
        scenario[config.parameter_name] = value
        scenarios.append(scenario)

    if self.verbose:
        print(f"Parameter sweep: {config.parameter_name}")
        print(f"  Values: {config.values}")

    return self.run_scenarios(scenarios, duration, initial_state, **kwargs)

run_scenarios(scenarios, duration, initial_state, dt=1.0 / 24.0, compute_metrics=True, save_time_series=False)

Run multiple scenarios in parallel.

Source code in pyadm1/simulation/parallel.py

def run_scenarios(
    self,
    scenarios: List[Dict[str, Any]],
    duration: float,
    initial_state: List[float],
    dt: float = 1.0 / 24.0,
    compute_metrics: bool = True,
    save_time_series: bool = False,
) -> List[ScenarioResult]:
    """Run multiple scenarios in parallel."""
    if self.verbose:
        print(f"Starting parallel simulation with {len(scenarios)} scenarios")
        print(f"Using {self.n_workers} worker processes")

    start_time = time.time()

    worker_func = partial(
        _run_single_scenario,
        adm1_config=self._serialize_adm1(),
        duration=duration,
        initial_state=initial_state,
        dt=dt,
        compute_metrics=compute_metrics,
        save_time_series=save_time_series,
        verbose=self.verbose,
    )

    scenarios_with_ids = [(i, scenario) for i, scenario in enumerate(scenarios)]

    use_sequential = self.n_workers <= 1 or len(scenarios_with_ids) <= 1

    if use_sequential:
        results = []
        for i, scenario_data in enumerate(scenarios_with_ids):
            results.append(worker_func(scenario_data))
            if self.verbose and ((i + 1) % 10 == 0 or (i + 1) == len(scenarios)):
                print(f"  Completed {i + 1}/{len(scenarios)} scenarios")
    else:
        ctx = _get_mp_context()
        with ctx.Pool(processes=self.n_workers) as pool:
            if self.verbose:
                results = []
                for i, result in enumerate(pool.imap(worker_func, scenarios_with_ids)):
                    results.append(result)
                    if (i + 1) % 10 == 0 or (i + 1) == len(scenarios):
                        print(f"  Completed {i + 1}/{len(scenarios)} scenarios")
            else:
                results = pool.map(worker_func, scenarios_with_ids)

    elapsed_time = time.time() - start_time

    if self.verbose:
        n_success = sum(1 for r in results if r.success)
        print("\nSimulation complete:")
        print(f"  Total scenarios: {len(scenarios)}")
        print(f"  Successful: {n_success}")
        print(f"  Failed: {len(scenarios) - n_success}")
        print(f"  Total time: {elapsed_time:.1f} seconds")
        print(f"  Average time per scenario: {elapsed_time / len(scenarios):.2f} seconds")

    return results

summarize_results(results, metrics=None)

Compute summary statistics across multiple scenarios.

Source code in pyadm1/simulation/parallel.py

def summarize_results(self, results: List[ScenarioResult], metrics: Optional[List[str]] = None) -> Dict[str, Any]:
    """Compute summary statistics across multiple scenarios."""
    successful = [r for r in results if r.success]
    n_scenarios = len(results)
    n_successful = len(successful)

    summary = {
        "n_scenarios": n_scenarios,
        "n_successful": n_successful,
        "n_failed": n_scenarios - n_successful,
        "success_rate": n_successful / n_scenarios if n_scenarios > 0 else 0,
        "metrics": {},
    }

    if not successful:
        if n_scenarios > 0:
            summary["error"] = "No successful scenarios"
        else:
            summary["error"] = "No scenarios to summarize"
        return summary

    if metrics is None:
        metrics = list(successful[0].metrics.keys())

    for metric_name in metrics:
        values = [r.metrics.get(metric_name, np.nan) for r in successful]
        values = [v for v in values if not np.isnan(v)]

        if values:
            summary["metrics"][metric_name] = {
                "mean": np.mean(values),
                "std": np.std(values),
                "min": np.min(values),
                "max": np.max(values),
                "median": np.median(values),
                "q25": np.percentile(values, 25),
                "q75": np.percentile(values, 75),
            }

    return summary

pyadm1.simulation.parallel.ScenarioResult dataclass

Result from a single simulation scenario.

Source code in pyadm1/simulation/parallel.py

@dataclass
class ScenarioResult:
    """Result from a single simulation scenario."""

    scenario_id: int
    parameters: Dict[str, Any]
    success: bool
    duration: float
    final_state: Optional[List[float]] = None
    time_series: Optional[Dict[str, List[float]]] = None
    metrics: Dict[str, float] = field(default_factory=dict)
    error: Optional[str] = None
    execution_time: float = 0.0

pyadm1.simulation.parallel.ParameterSweepConfig dataclass

Configuration for a single-parameter sweep.

Source code in pyadm1/simulation/parallel.py

@dataclass
class ParameterSweepConfig:
    """Configuration for a single-parameter sweep."""

    parameter_name: str
    values: List[float]
    other_params: Dict[str, Any] = field(default_factory=dict)

pyadm1.simulation.parallel.MonteCarloConfig dataclass

Configuration for Monte Carlo simulation.

Source code in pyadm1/simulation/parallel.py

@dataclass
class MonteCarloConfig:
    """Configuration for Monte Carlo simulation."""

    n_samples: int
    parameter_distributions: Dict[str, Tuple[float, float]]
    fixed_params: Dict[str, Any] = field(default_factory=dict)
    seed: Optional[int] = None