Skip to content

Simulation API

pyadm1.simulation.simulator.Simulator

Handles ADM1 simulation runs with various configurations.

This class provides high-level interfaces for running ADM1 simulations, including single runs and multi-scenario optimization for substrate feed determination.

Attributes:

Name Type Description
adm1 ADM1

ADM1 model instance
solver ODESolver

ODE solver instance
Example

simulator = Simulator(adm1) result = simulator.simulate_AD_plant([0, 10], initial_state)

Source code in pyadm1/simulation/simulator.py 
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
class Simulator:
    """
    Handles ADM1 simulation runs with various configurations.

    This class provides high-level interfaces for running ADM1 simulations,
    including single runs and multi-scenario optimization for substrate feed
    determination.

    Attributes:
        adm1: ADM1 model instance
        solver: ODE solver instance

    Example:
        >>> simulator = Simulator(adm1)
        >>> result = simulator.simulate_AD_plant([0, 10], initial_state)
    """

    def __init__(self, adm1: ADM1, solver: ODESolver = None) -> None:
        """
        Initialize simulator with ADM1 model instance.

        Args:
            adm1: ADM1 model instance
            solver: Optional ODE solver. If None, creates default BDF solver
        """
        self._adm1 = adm1
        self._solver = solver or create_solver(method="BDF", rtol=1e-6, atol=1e-8)

    def determine_best_feed_by_n_sims(
        self,
        state_zero: List[float],
        Q: List[float],
        Qch4sp: float,
        feeding_freq: int,
        n: int = 13,
    ) -> Tuple[float, float, List[float], float, float, List[float], float, float, float, float]:
        """
        Determine optimal substrate feed by running n simulations.

        Runs n simulations with varying substrate feed rates around Q and
        returns the feed rate yielding methane production closest to setpoint.

        The first simulation uses Q, the 2nd and 3rd use Q ± 1.5 m³/d,
        and remaining simulations use random variations.

        Args:
            state_zero: Initial ADM1 state vector (37 elements)
            Q: Initial volumetric flow rates [m³/d], e.g. [15, 10, 0, ...]
            Qch4sp: Methane flow rate setpoint [m³/d]
            feeding_freq: Feeding frequency [hours]
            n: Number of simulations to run (default: 13, minimum: 3)

        Returns:
            Tuple containing:
                - Q_Gas_7d_best: Best biogas production after 7 days [m³/d]
                - Q_CH4_7d_best: Best methane production after 7 days [m³/d]
                - Qbest: Best substrate feed rates [m³/d]
                - Q_Gas_7d_initial: Initial biogas production after 7 days [m³/d]
                - Q_CH4_7d_initial: Initial methane production after 7 days [m³/d]
                - Q_initial: Initial substrate feed rates [m³/d]
                - q_gas_best_2d: Best biogas after feeding_freq/24 days [m³/d]
                - q_ch4_best_2d: Best methane after feeding_freq/24 days [m³/d]
                - q_gas_2d: Initial biogas after feeding_freq/24 days [m³/d]
                - q_ch4_2d: Initial methane after feeding_freq/24 days [m³/d]

        Example:
            >>> result = simulator.determine_best_feed_by_n_sims(
            ...     state, [15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 900, 48, n=13
            ... )
            >>> Q_best = result[2]
        """
        if n < 3:
            raise ValueError("n must be at least 3")

        Q_CH4_7d = [0.0] * n
        Q_Gas_7d = [0.0] * n

        # Generate substrate feed mixtures
        Qnew = self._adm1.feedstock.get_substrate_feed_mixtures(Q, n)

        # Run n simulations
        for i, q in enumerate(Qnew):
            Q_Gas_7d[i], Q_CH4_7d[i] = self._simulate_without_saving_state([0, 7], state_zero, q)

        # Find scenario closest to methane setpoint
        ii = np.argmin([np.sum((q - Qch4sp) ** 2) for q in Q_CH4_7d])
        Qbest = Qnew[ii]

        # Simulate with initial Q for feeding_freq/24 days
        q_gas_2d, q_ch4_2d = self._simulate_without_saving_state([0, feeding_freq / 24], state_zero, Qnew[0])

        # Simulate with best Q for feeding_freq/24 days
        q_gas_best_2d, q_ch4_best_2d = self._simulate_without_saving_state([0, feeding_freq / 24], state_zero, Qbest)

        return (
            Q_Gas_7d[ii][-1] if hasattr(Q_Gas_7d[ii], "__iter__") else Q_Gas_7d[ii],
            Q_CH4_7d[ii][-1] if hasattr(Q_CH4_7d[ii], "__iter__") else Q_CH4_7d[ii],
            Qbest,
            Q_Gas_7d[0][-1] if hasattr(Q_Gas_7d[0], "__iter__") else Q_Gas_7d[0],
            Q_CH4_7d[0][-1] if hasattr(Q_CH4_7d[0], "__iter__") else Q_CH4_7d[0],
            Qnew[0],
            q_gas_best_2d[-1] if hasattr(q_gas_best_2d, "__iter__") else q_gas_best_2d,
            q_ch4_best_2d[-1] if hasattr(q_ch4_best_2d, "__iter__") else q_ch4_best_2d,
            q_gas_2d[-1] if hasattr(q_gas_2d, "__iter__") else q_gas_2d,
            q_ch4_2d[-1] if hasattr(q_ch4_2d, "__iter__") else q_ch4_2d,
        )

    def simulate_AD_plant(self, tstep: List[float], state_zero: List[float]) -> List[float]:
        """
        Simulate ADM1 for specified time span and return final state.

        This is the main simulation method that integrates the ADM1 ODEs
        and tracks process values for operator information.

        Args:
            tstep: Time span [t_start, t_end] in days
            state_zero: Initial ADM1 state vector (37 elements)

        Returns:
            Final ADM1 state vector after simulation (37 elements)

        Example:
            >>> final_state = simulator.simulate_AD_plant([0, 1], initial_state)
            >>> print(f"Final pH: {final_state[...])
        """
        final_state = self._simulate_and_return_final_state(tstep, state_zero)

        # Print process parameters for monitoring
        self._adm1.print_params_at_current_state(final_state)

        return final_state

    def _simulate_without_saving_state(
        self, tstep: List[float], state_zero: List[float], Q: List[float]
    ) -> Tuple[float, float]:
        """
        Run simulation without saving final state.

        Used internally for optimization scenarios where we only need
        gas production values.

        Args:
            tstep: Time span [t_start, t_end] in days
            state_zero: Initial ADM1 state vector
            Q: Volumetric flow rates [m³/d]

        Returns:
            Tuple of (q_gas, q_ch4) - biogas and methane production rates [m³/d]
        """
        # Create ADM1 input stream
        self._adm1.create_influent(Q, 0)

        # Integrate ODEs
        result = self._solver.solve(fun=self._adm1.ADM1_ODE, t_span=tstep, y0=state_zero)

        # Extract final gas phase state
        final_state = result.y[:, -1]
        pi_Sh2, pi_Sch4, pi_Sco2, pTOTAL = final_state[33:37]

        # Calculate gas production rates
        q_gas, q_ch4, q_co2, p_gas = self._adm1.calc_gas(pi_Sh2, pi_Sch4, pi_Sco2, pTOTAL)

        return q_gas, q_ch4

    def _simulate_and_return_final_state(self, tstep: List[float], state_zero: List[float]) -> List[float]:
        """
        Simulate and return final state vector.

        Args:
            tstep: Time span [t_start, t_end] in days
            state_zero: Initial ADM1 state vector

        Returns:
            Final ADM1 state vector (37 elements)
        """
        # Integrate ODEs using solver
        result = self._solver.solve(fun=self._adm1.ADM1_ODE, t_span=tstep, y0=state_zero)

        # Extract final state
        final_state = result.y[:, -1].tolist()

        return final_state

    @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

__init__(adm1, solver=None)

Initialize simulator with ADM1 model instance.

Parameters:

Name Type Description Default
adm1 ADM1

ADM1 model instance

 required
solver ODESolver

Optional ODE solver. If None, creates default BDF solver

 None
Source code in pyadm1/simulation/simulator.py 
46
47
48
49
50
51
52
53
54
55
def __init__(self, adm1: ADM1, solver: ODESolver = None) -> None:
    """
    Initialize simulator with ADM1 model instance.

    Args:
        adm1: ADM1 model instance
        solver: Optional ODE solver. If None, creates default BDF solver
    """
    self._adm1 = adm1
    self._solver = solver or create_solver(method="BDF", rtol=1e-6, atol=1e-8)

determine_best_feed_by_n_sims(state_zero, Q, Qch4sp, feeding_freq, n=13)

Determine optimal substrate feed by running n simulations.

Runs n simulations with varying substrate feed rates around Q and returns the feed rate yielding methane production closest to setpoint.

The first simulation uses Q, the 2nd and 3rd use Q ± 1.5 m³/d, and remaining simulations use random variations.

Parameters:

Name Type Description Default
state_zero List[float]

Initial ADM1 state vector (37 elements)

 required
Q List[float]

Initial volumetric flow rates [m³/d], e.g. [15, 10, 0, ...]

 required
Qch4sp float

Methane flow rate setpoint [m³/d]

 required
feeding_freq int

Feeding frequency [hours]

 required
n int

Number of simulations to run (default: 13, minimum: 3)

 13

Returns:

Type Description
Tuple[float, float, List[float], float, float, List[float], float, float, float, float]

Tuple containing: - Q_Gas_7d_best: Best biogas production after 7 days [m³/d] - Q_CH4_7d_best: Best methane production after 7 days [m³/d] - Qbest: Best substrate feed rates [m³/d] - Q_Gas_7d_initial: Initial biogas production after 7 days [m³/d] - Q_CH4_7d_initial: Initial methane production after 7 days [m³/d] - Q_initial: Initial substrate feed rates [m³/d] - q_gas_best_2d: Best biogas after feeding_freq/24 days [m³/d] - q_ch4_best_2d: Best methane after feeding_freq/24 days [m³/d] - q_gas_2d: Initial biogas after feeding_freq/24 days [m³/d] - q_ch4_2d: Initial methane after feeding_freq/24 days [m³/d]
Example

result = simulator.determine_best_feed_by_n_sims( ... state, [15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 900, 48, n=13 ... ) Q_best = result[2]

Source code in pyadm1/simulation/simulator.py 
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def determine_best_feed_by_n_sims(
    self,
    state_zero: List[float],
    Q: List[float],
    Qch4sp: float,
    feeding_freq: int,
    n: int = 13,
) -> Tuple[float, float, List[float], float, float, List[float], float, float, float, float]:
    """
    Determine optimal substrate feed by running n simulations.

    Runs n simulations with varying substrate feed rates around Q and
    returns the feed rate yielding methane production closest to setpoint.

    The first simulation uses Q, the 2nd and 3rd use Q ± 1.5 m³/d,
    and remaining simulations use random variations.

    Args:
        state_zero: Initial ADM1 state vector (37 elements)
        Q: Initial volumetric flow rates [m³/d], e.g. [15, 10, 0, ...]
        Qch4sp: Methane flow rate setpoint [m³/d]
        feeding_freq: Feeding frequency [hours]
        n: Number of simulations to run (default: 13, minimum: 3)

    Returns:
        Tuple containing:
            - Q_Gas_7d_best: Best biogas production after 7 days [m³/d]
            - Q_CH4_7d_best: Best methane production after 7 days [m³/d]
            - Qbest: Best substrate feed rates [m³/d]
            - Q_Gas_7d_initial: Initial biogas production after 7 days [m³/d]
            - Q_CH4_7d_initial: Initial methane production after 7 days [m³/d]
            - Q_initial: Initial substrate feed rates [m³/d]
            - q_gas_best_2d: Best biogas after feeding_freq/24 days [m³/d]
            - q_ch4_best_2d: Best methane after feeding_freq/24 days [m³/d]
            - q_gas_2d: Initial biogas after feeding_freq/24 days [m³/d]
            - q_ch4_2d: Initial methane after feeding_freq/24 days [m³/d]

    Example:
        >>> result = simulator.determine_best_feed_by_n_sims(
        ...     state, [15, 10, 0, 0, 0, 0, 0, 0, 0, 0], 900, 48, n=13
        ... )
        >>> Q_best = result[2]
    """
    if n < 3:
        raise ValueError("n must be at least 3")

    Q_CH4_7d = [0.0] * n
    Q_Gas_7d = [0.0] * n

    # Generate substrate feed mixtures
    Qnew = self._adm1.feedstock.get_substrate_feed_mixtures(Q, n)

    # Run n simulations
    for i, q in enumerate(Qnew):
        Q_Gas_7d[i], Q_CH4_7d[i] = self._simulate_without_saving_state([0, 7], state_zero, q)

    # Find scenario closest to methane setpoint
    ii = np.argmin([np.sum((q - Qch4sp) ** 2) for q in Q_CH4_7d])
    Qbest = Qnew[ii]

    # Simulate with initial Q for feeding_freq/24 days
    q_gas_2d, q_ch4_2d = self._simulate_without_saving_state([0, feeding_freq / 24], state_zero, Qnew[0])

    # Simulate with best Q for feeding_freq/24 days
    q_gas_best_2d, q_ch4_best_2d = self._simulate_without_saving_state([0, feeding_freq / 24], state_zero, Qbest)

    return (
        Q_Gas_7d[ii][-1] if hasattr(Q_Gas_7d[ii], "__iter__") else Q_Gas_7d[ii],
        Q_CH4_7d[ii][-1] if hasattr(Q_CH4_7d[ii], "__iter__") else Q_CH4_7d[ii],
        Qbest,
        Q_Gas_7d[0][-1] if hasattr(Q_Gas_7d[0], "__iter__") else Q_Gas_7d[0],
        Q_CH4_7d[0][-1] if hasattr(Q_CH4_7d[0], "__iter__") else Q_CH4_7d[0],
        Qnew[0],
        q_gas_best_2d[-1] if hasattr(q_gas_best_2d, "__iter__") else q_gas_best_2d,
        q_ch4_best_2d[-1] if hasattr(q_ch4_best_2d, "__iter__") else q_ch4_best_2d,
        q_gas_2d[-1] if hasattr(q_gas_2d, "__iter__") else q_gas_2d,
        q_ch4_2d[-1] if hasattr(q_ch4_2d, "__iter__") else q_ch4_2d,
    )

simulate_AD_plant(tstep, state_zero)

Simulate ADM1 for specified time span and return final state.

This is the main simulation method that integrates the ADM1 ODEs and tracks process values for operator information.

Parameters:

Name Type Description Default
tstep List[float]

Time span [t_start, t_end] in days

 required
state_zero List[float]

Initial ADM1 state vector (37 elements)

 required

Returns:

Type Description
List[float]

Final ADM1 state vector after simulation (37 elements)
Example

final_state = simulator.simulate_AD_plant([0, 1], initial_state) print(f"Final pH: {final_state[...])

Source code in pyadm1/simulation/simulator.py 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def simulate_AD_plant(self, tstep: List[float], state_zero: List[float]) -> List[float]:
    """
    Simulate ADM1 for specified time span and return final state.

    This is the main simulation method that integrates the ADM1 ODEs
    and tracks process values for operator information.

    Args:
        tstep: Time span [t_start, t_end] in days
        state_zero: Initial ADM1 state vector (37 elements)

    Returns:
        Final ADM1 state vector after simulation (37 elements)

    Example:
        >>> final_state = simulator.simulate_AD_plant([0, 1], initial_state)
        >>> print(f"Final pH: {final_state[...])
    """
    final_state = self._simulate_and_return_final_state(tstep, state_zero)

    # Print process parameters for monitoring
    self._adm1.print_params_at_current_state(final_state)

    return final_state

pyadm1.simulation.parallel.ParallelSimulator

Parallel simulator for running multiple ADM1 scenarios concurrently.

Uses multiprocessing to distribute scenarios across CPU cores for efficient parameter sweeps, sensitivity analysis, and Monte Carlo simulations.

Attributes:

Name Type Description
adm1

Base ADM1 model instance (will be copied for each worker)
n_workers

Number of parallel worker processes
verbose

Enable progress reporting
Example

parallel = ParallelSimulator(adm1, n_workers=4, verbose=True) results = parallel.run_scenarios(scenarios, duration=30)

Source code in pyadm1/simulation/parallel.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
class ParallelSimulator:
    """
    Parallel simulator for running multiple ADM1 scenarios concurrently.

    Uses multiprocessing to distribute scenarios across CPU cores for efficient
    parameter sweeps, sensitivity analysis, and Monte Carlo simulations.

    Attributes:
        adm1: Base ADM1 model instance (will be copied for each worker)
        n_workers: Number of parallel worker processes
        verbose: Enable progress reporting

    Example:
        >>> parallel = ParallelSimulator(adm1, n_workers=4, verbose=True)
        >>> results = parallel.run_scenarios(scenarios, duration=30)
    """

    def __init__(self, adm1: "ADM1", n_workers: Optional[int] = None, verbose: bool = True):
        """
        Initialize parallel simulator.

        Args:
            adm1: ADM1 model instance
            n_workers: Number of worker processes (default: CPU count - 1)
            verbose: Enable progress output
        """
        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 simulation scenarios in parallel.

        Each scenario is a dictionary containing parameter values and substrate
        feed rates. The simulator will run all scenarios concurrently and
        collect results.

        Args:
            scenarios: List of scenario dictionaries with parameters
            duration: Simulation duration [days]
            initial_state: Initial ADM1 state vector
            dt: Time step [days]
            compute_metrics: Calculate performance metrics
            save_time_series: Save full time series data

        Returns:
            List of ScenarioResult objects

        Example:
            >>> scenarios = [
            ...     {"k_dis": 0.5, "Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]},
            ...     {"k_dis": 0.6, "Q": [20, 10, 0, 0, 0, 0, 0, 0, 0, 0]},
            ... ]
            >>> results = parallel.run_scenarios(scenarios, duration=30, initial_state=state)
        """
        if self.verbose:
            print(f"Starting parallel simulation with {len(scenarios)} scenarios")
            print(f"Using {self.n_workers} worker processes")

        start_time = time.time()

        # Create worker function with fixed parameters
        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,
        )

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

        # Avoid spawning a process pool for trivial workloads. This also prevents
        # pythonnet/.NET deadlocks seen in some CI environments with forked workers.
        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()
            # Run scenarios in parallel
            # Keep workers alive across multiple scenarios to avoid repeated
            # python startup / pythonnet initialization overhead.
            with ctx.Pool(processes=self.n_workers) as pool:
                if self.verbose:
                    # Use imap for progress tracking
                    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 parameter sweep for a single parameter.

        Tests multiple values of one parameter while keeping others fixed.

        Args:
            config: ParameterSweepConfig with parameter and values
            duration: Simulation duration [days]
            initial_state: Initial ADM1 state vector
            **kwargs: Additional arguments for run_scenarios

        Returns:
            List of ScenarioResult objects

        Example:
            >>> config = ParameterSweepConfig(
            ...     parameter_name="k_dis",
            ...     values=[0.3, 0.4, 0.5, 0.6, 0.7],
            ...     other_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]}
            ... )
            >>> results = parallel.parameter_sweep(config, duration=30, initial_state=state)
        """
        # Generate scenarios
        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 multi-parameter sweep (full factorial design).

        Tests all combinations of provided parameter values.

        Args:
            parameter_configs: Dict mapping parameter names to value lists
            duration: Simulation duration [days]
            initial_state: Initial ADM1 state vector
            fixed_params: Parameters to keep fixed
            **kwargs: Additional arguments for run_scenarios

        Returns:
            List of ScenarioResult objects

        Example:
            >>> parameter_configs = {
            ...     "k_dis": [0.4, 0.5, 0.6],
            ...     "Y_su": [0.09, 0.10, 0.11]
            ... }
            >>> results = parallel.multi_parameter_sweep(
            ...     parameter_configs,
            ...     duration=30,
            ...     initial_state=state,
            ...     fixed_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]}
            ... )
        """
        fixed_params = fixed_params or {}

        # Generate all combinations using recursive approach
        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.

        Samples parameters from normal distributions and runs multiple scenarios
        to quantify uncertainty in predictions.

        Args:
            config: MonteCarloConfig with distributions and sample count
            duration: Simulation duration [days]
            initial_state: Initial ADM1 state vector
            **kwargs: Additional arguments for run_scenarios

        Returns:
            List of ScenarioResult objects

        Example:
            >>> config = MonteCarloConfig(
            ...     n_samples=100,
            ...     parameter_distributions={
            ...         "k_dis": (0.5, 0.05),  # mean=0.5, std=0.05
            ...         "Y_su": (0.10, 0.01)
            ...     },
            ...     fixed_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]},
            ...     seed=42
            ... )
            >>> results = parallel.monte_carlo(config, duration=30, initial_state=state)
        """
        # Set random seed for reproducibility
        if config.seed is not None:
            np.random.seed(config.seed)

        # Generate scenarios
        scenarios = []
        for i in range(config.n_samples):
            scenario = config.fixed_params.copy()

            # Sample each parameter from its distribution
            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]:
        """
        Summarize results from multiple scenarios.

        Computes summary statistics for each metric across all successful scenarios.

        Args:
            results: List of ScenarioResult objects
            metrics: List of metric names to summarize (default: all)

        Returns:
            Dictionary with summary statistics

        Example:
            ```python
            summary = parallel.summarize_results(results)
            ch4_stats = summary['metrics']['Q_ch4']
            print(f"Mean CH4: {ch4_stats['avg']:.1f} m³/d")
            ```
        """
        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

        # Determine which metrics to summarize
        if metrics is None:
            # Use all metrics from first result
            metrics = list(successful[0].metrics.keys())

        # Compute statistics for each metric
        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 for passing to worker processes.

        Returns:
            Dictionary with ADM1 configuration
        """
        return {
            "V_liq": self.adm1.V_liq,
            "V_gas": self.adm1._V_gas,
            "T_ad": self.adm1._T_ad,
            "feedstock_config": {
                "feeding_freq": 48,  # Default from feedstock
            },
        }

    @staticmethod
    def _generate_combinations(value_lists: List[List[Any]]) -> List[Tuple]:
        """
        Generate all combinations from lists of values (Cartesian product).

        Args:
            value_lists: List of lists containing values

        Returns:
            List of tuples with all combinations
        """
        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)

Initialize parallel simulator.

Parameters:

Name Type Description Default
adm1 ADM1

ADM1 model instance

 required
n_workers Optional[int]

Number of worker processes (default: CPU count - 1)

 None
verbose bool

Enable progress output

 True
Source code in pyadm1/simulation/parallel.py 
177
178
179
180
181
182
183
184
185
186
187
188
def __init__(self, adm1: "ADM1", n_workers: Optional[int] = None, verbose: bool = True):
    """
    Initialize parallel simulator.

    Args:
        adm1: ADM1 model instance
        n_workers: Number of worker processes (default: CPU count - 1)
        verbose: Enable progress output
    """
    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.

Samples parameters from normal distributions and runs multiple scenarios to quantify uncertainty in predictions.

Parameters:

Name Type Description Default
config MonteCarloConfig

MonteCarloConfig with distributions and sample count

 required
duration float

Simulation duration [days]

 required
initial_state List[float]

Initial ADM1 state vector

 required
**kwargs Any

Additional arguments for run_scenarios

 {}

Returns:

Type Description
List[ScenarioResult]

List of ScenarioResult objects
Example

config = MonteCarloConfig( ... n_samples=100, ... parameter_distributions={ ... "k_dis": (0.5, 0.05), # mean=0.5, std=0.05 ... "Y_su": (0.10, 0.01) ... }, ... fixed_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]}, ... seed=42 ... ) results = parallel.monte_carlo(config, duration=30, initial_state=state)

Source code in pyadm1/simulation/parallel.py 
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
def monte_carlo(
    self,
    config: MonteCarloConfig,
    duration: float,
    initial_state: List[float],
    **kwargs: Any,
) -> List[ScenarioResult]:
    """
    Run Monte Carlo simulation with parameter uncertainty.

    Samples parameters from normal distributions and runs multiple scenarios
    to quantify uncertainty in predictions.

    Args:
        config: MonteCarloConfig with distributions and sample count
        duration: Simulation duration [days]
        initial_state: Initial ADM1 state vector
        **kwargs: Additional arguments for run_scenarios

    Returns:
        List of ScenarioResult objects

    Example:
        >>> config = MonteCarloConfig(
        ...     n_samples=100,
        ...     parameter_distributions={
        ...         "k_dis": (0.5, 0.05),  # mean=0.5, std=0.05
        ...         "Y_su": (0.10, 0.01)
        ...     },
        ...     fixed_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]},
        ...     seed=42
        ... )
        >>> results = parallel.monte_carlo(config, duration=30, initial_state=state)
    """
    # Set random seed for reproducibility
    if config.seed is not None:
        np.random.seed(config.seed)

    # Generate scenarios
    scenarios = []
    for i in range(config.n_samples):
        scenario = config.fixed_params.copy()

        # Sample each parameter from its distribution
        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 multi-parameter sweep (full factorial design).

Tests all combinations of provided parameter values.

Parameters:

Name Type Description Default
parameter_configs Dict[str, List[float]]

Dict mapping parameter names to value lists

 required
duration float

Simulation duration [days]

 required
initial_state List[float]

Initial ADM1 state vector

 required
fixed_params Optional[Dict[str, Any]]

Parameters to keep fixed

 None
**kwargs Any

Additional arguments for run_scenarios

 {}

Returns:

Type Description
List[ScenarioResult]

List of ScenarioResult objects
Example

parameter_configs = { ... "k_dis": [0.4, 0.5, 0.6], ... "Y_su": [0.09, 0.10, 0.11] ... } results = parallel.multi_parameter_sweep( ... parameter_configs, ... duration=30, ... initial_state=state, ... fixed_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]} ... )

Source code in pyadm1/simulation/parallel.py 
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
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 multi-parameter sweep (full factorial design).

    Tests all combinations of provided parameter values.

    Args:
        parameter_configs: Dict mapping parameter names to value lists
        duration: Simulation duration [days]
        initial_state: Initial ADM1 state vector
        fixed_params: Parameters to keep fixed
        **kwargs: Additional arguments for run_scenarios

    Returns:
        List of ScenarioResult objects

    Example:
        >>> parameter_configs = {
        ...     "k_dis": [0.4, 0.5, 0.6],
        ...     "Y_su": [0.09, 0.10, 0.11]
        ... }
        >>> results = parallel.multi_parameter_sweep(
        ...     parameter_configs,
        ...     duration=30,
        ...     initial_state=state,
        ...     fixed_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]}
        ... )
    """
    fixed_params = fixed_params or {}

    # Generate all combinations using recursive approach
    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 parameter sweep for a single parameter.

Tests multiple values of one parameter while keeping others fixed.

Parameters:

Name Type Description Default
config ParameterSweepConfig

ParameterSweepConfig with parameter and values

 required
duration float

Simulation duration [days]

 required
initial_state List[float]

Initial ADM1 state vector

 required
**kwargs Any

Additional arguments for run_scenarios

 {}

Returns:

Type Description
List[ScenarioResult]

List of ScenarioResult objects
Example

config = ParameterSweepConfig( ... parameter_name="k_dis", ... values=[0.3, 0.4, 0.5, 0.6, 0.7], ... other_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]} ... ) results = parallel.parameter_sweep(config, duration=30, initial_state=state)

Source code in pyadm1/simulation/parallel.py 
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
def parameter_sweep(
    self,
    config: ParameterSweepConfig,
    duration: float,
    initial_state: List[float],
    **kwargs: Any,
) -> List[ScenarioResult]:
    """
    Run parameter sweep for a single parameter.

    Tests multiple values of one parameter while keeping others fixed.

    Args:
        config: ParameterSweepConfig with parameter and values
        duration: Simulation duration [days]
        initial_state: Initial ADM1 state vector
        **kwargs: Additional arguments for run_scenarios

    Returns:
        List of ScenarioResult objects

    Example:
        >>> config = ParameterSweepConfig(
        ...     parameter_name="k_dis",
        ...     values=[0.3, 0.4, 0.5, 0.6, 0.7],
        ...     other_params={"Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]}
        ... )
        >>> results = parallel.parameter_sweep(config, duration=30, initial_state=state)
    """
    # Generate scenarios
    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 simulation scenarios in parallel.

Each scenario is a dictionary containing parameter values and substrate feed rates. The simulator will run all scenarios concurrently and collect results.

Parameters:

Name Type Description Default
scenarios List[Dict[str, Any]]

List of scenario dictionaries with parameters

 required
duration float

Simulation duration [days]

 required
initial_state List[float]

Initial ADM1 state vector

 required
dt float

Time step [days]

 1.0 / 24.0
compute_metrics bool

Calculate performance metrics

 True
save_time_series bool

Save full time series data

 False

Returns:

Type Description
List[ScenarioResult]

List of ScenarioResult objects
Example

scenarios = [ ... {"k_dis": 0.5, "Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]}, ... {"k_dis": 0.6, "Q": [20, 10, 0, 0, 0, 0, 0, 0, 0, 0]}, ... ] results = parallel.run_scenarios(scenarios, duration=30, initial_state=state)

Source code in pyadm1/simulation/parallel.py 
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
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 simulation scenarios in parallel.

    Each scenario is a dictionary containing parameter values and substrate
    feed rates. The simulator will run all scenarios concurrently and
    collect results.

    Args:
        scenarios: List of scenario dictionaries with parameters
        duration: Simulation duration [days]
        initial_state: Initial ADM1 state vector
        dt: Time step [days]
        compute_metrics: Calculate performance metrics
        save_time_series: Save full time series data

    Returns:
        List of ScenarioResult objects

    Example:
        >>> scenarios = [
        ...     {"k_dis": 0.5, "Q": [15, 10, 0, 0, 0, 0, 0, 0, 0, 0]},
        ...     {"k_dis": 0.6, "Q": [20, 10, 0, 0, 0, 0, 0, 0, 0, 0]},
        ... ]
        >>> results = parallel.run_scenarios(scenarios, duration=30, initial_state=state)
    """
    if self.verbose:
        print(f"Starting parallel simulation with {len(scenarios)} scenarios")
        print(f"Using {self.n_workers} worker processes")

    start_time = time.time()

    # Create worker function with fixed parameters
    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,
    )

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

    # Avoid spawning a process pool for trivial workloads. This also prevents
    # pythonnet/.NET deadlocks seen in some CI environments with forked workers.
    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()
        # Run scenarios in parallel
        # Keep workers alive across multiple scenarios to avoid repeated
        # python startup / pythonnet initialization overhead.
        with ctx.Pool(processes=self.n_workers) as pool:
            if self.verbose:
                # Use imap for progress tracking
                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)

Summarize results from multiple scenarios.

Computes summary statistics for each metric across all successful scenarios.

Parameters:

Name Type Description Default
results List[ScenarioResult]

List of ScenarioResult objects

 required
metrics Optional[List[str]]

List of metric names to summarize (default: all)

 None

Returns:

Type Description
Dict[str, Any]

Dictionary with summary statistics
Example 
summary = parallel.summarize_results(results)
ch4_stats = summary['metrics']['Q_ch4']
print(f"Mean CH4: {ch4_stats['avg']:.1f} m³/d")
Source code in pyadm1/simulation/parallel.py 
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
def summarize_results(self, results: List[ScenarioResult], metrics: Optional[List[str]] = None) -> Dict[str, Any]:
    """
    Summarize results from multiple scenarios.

    Computes summary statistics for each metric across all successful scenarios.

    Args:
        results: List of ScenarioResult objects
        metrics: List of metric names to summarize (default: all)

    Returns:
        Dictionary with summary statistics

    Example:
        ```python
        summary = parallel.summarize_results(results)
        ch4_stats = summary['metrics']['Q_ch4']
        print(f"Mean CH4: {ch4_stats['avg']:.1f} m³/d")
        ```
    """
    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

    # Determine which metrics to summarize
    if metrics is None:
        # Use all metrics from first result
        metrics = list(successful[0].metrics.keys())

    # Compute statistics for each metric
    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