Batch / Parallel Design of Experiments

Batch / Parallel Design of Experiments#

When experimental capacity is parallel (several reactors, wells, or robots that run together), designing experiments one at a time wastes information. This tutorial demonstrates discopt.doe.batch_optimal_experiment, which selects \(N\) experimental conditions jointly under D/A/E/ME optimality.

Mathematical background. For independent experiments the Fisher information matrix is additive,

\[ \text{FIM}_{\text{total}}(d_1, \dots, d_N) = \sum_{i=1}^N J(\theta, d_i)^\top \Sigma^{-1} J(\theta, d_i) + \text{FIM}_{\text{prior}}, \]

and the joint D-optimal batch maximises \(\log\det \text{FIM}_{\text{total}}\) over the stacked design vector \((d_1, \dots, d_N)\). This formulation follows Galvanin et al. [2007]; we also support a cheaper greedy strategy (pick one, fold into the prior, repeat) and a distance-penalized variant. For robust / expected-information variants the reader is referred to Sandrin et al. [2025]. For an alternative route via parallel Bayesian optimization rather than FIM-based MBDoE, see González and Zavala [2023].

import os

os.environ.setdefault("JAX_PLATFORMS", "cpu")
os.environ.setdefault("JAX_ENABLE_X64", "1")

import discopt.modeling as dm
import numpy as np
from discopt.doe import (
    BatchStrategy,
    DesignCriterion,
    batch_optimal_experiment,
    sequential_doe,
)
from discopt.estimate import Experiment, ExperimentModel

Toy experiment#

A first-order decay \(y = A \exp(-k t)\) with one design variable (sampling time \(t\)) and one unknown parameter (\(k\)).

class ExpDecay(Experiment):
    def create_model(self, **kwargs):
        m = dm.Model("expdecay")
        k = m.continuous("k", lb=0.01, ub=5)
        t = m.continuous("t", lb=0.1, ub=10)
        return ExperimentModel(
            model=m,
            unknown_parameters={"k": k},
            design_inputs={"t": t},
            responses={"y": 5.0 * dm.exp(-k * t)},
            measurement_error={"y": 0.05},
        )


exp = ExpDecay()
param_values = {"k": 0.5}
bounds = {"t": (0.1, 10.0)}

Greedy batch#

Pick one design, fold its FIM into the running prior, repeat. Cheap and incremental, and it reuses the single-design optimal_experiment search.

greedy = batch_optimal_experiment(
    exp,
    param_values=param_values,
    design_bounds=bounds,
    n_experiments=4,
    strategy=BatchStrategy.GREEDY,
    criterion=DesignCriterion.D_OPTIMAL,
)
print(greedy.summary())
print()
print("per-round criterion:", greedy.per_round_criterion)

Batch Optimal Design (N=4, strategy=greedy)
============================================================
  Experiment 1:
                  t = 2
  Experiment 2:
                  t = 2
  Experiment 3:
                  t = 1.99996
  Experiment 4:
                  t = 2.00002

  D-opt (log det joint FIM) = 9.983
  A-opt (trace joint FIM⁻¹) = 4.618e-05
  E-opt (min eigenval)      = 2.165e+04
  Condition number          = 1

  SE(k) = 0.006796

per-round criterion: [8.596634833096028, 9.289782013655982, 9.69524712162411, 9.982929194076636]

Joint batch (Galvanina-style)#

Stack the \(N\) design vectors and optimize them simultaneously with scipy L-BFGS-B over the joint criterion. Slower than greedy on nonlinear experiments but equal or better in information content. The formulation follows Galvanin et al. [2007].

joint = batch_optimal_experiment(
    exp,
    param_values=param_values,
    design_bounds=bounds,
    n_experiments=4,
    strategy=BatchStrategy.JOINT,
    criterion=DesignCriterion.D_OPTIMAL,
)
print(joint.summary())
print()
print(f"greedy log det FIM = {greedy.criterion_value:.4g}")
print(f"joint  log det FIM = {joint.criterion_value:.4g}")

Batch Optimal Design (N=4, strategy=joint)
============================================================
  Experiment 1:
                  t = 1.99995
  Experiment 2:
                  t = 1.99997
  Experiment 3:
                  t = 2.00001
  Experiment 4:
                  t = 2

  D-opt (log det joint FIM) = 9.983
  A-opt (trace joint FIM⁻¹) = 4.618e-05
  E-opt (min eigenval)      = 2.165e+04
  Condition number          = 1

  SE(k) = 0.006796

greedy log det FIM = 9.983
joint  log det FIM = 9.983

Penalized batch#

Greedy with a normalized min-distance filter: useful when you want spatial diversity across the batch (e.g. to avoid clustering probes even if the criterion is locally indifferent).

penalized = batch_optimal_experiment(
    exp,
    param_values=param_values,
    design_bounds=bounds,
    n_experiments=4,
    strategy=BatchStrategy.PENALIZED,
    criterion=DesignCriterion.D_OPTIMAL,
    min_distance=0.15,
)
print(penalized.summary())

Batch Optimal Design (N=4, strategy=penalized)
============================================================
  Experiment 1:
                  t = 2
  Experiment 2:
                  t = 4.25778
  Experiment 3:
                  t = 6.20369
  Experiment 4:
                  t = 7.81546

  D-opt (log det joint FIM) = 9.105
  A-opt (trace joint FIM⁻¹) = 0.0001111
  E-opt (min eigenval)      = 9004
  Condition number          = 1

  SE(k) = 0.01054

Sequential batch MBDoE#

sequential_doe(experiments_per_round=N, batch_strategy=...) alternates parameter estimation and batch design, which is the pattern typically deployed in an automated laboratory. Each round designs \(N\) experiments jointly and the runner callback is invoked once per batch member.

k_true = 0.7


def runner(design):
    t = design["t"]
    return {"y": float(5.0 * np.exp(-k_true * t))}


history = sequential_doe(
    exp,
    initial_data={"y": np.array([5.0, 3.0])},
    initial_guess={"k": 0.5},
    design_bounds=bounds,
    n_rounds=3,
    experiments_per_round=2,
    batch_strategy=BatchStrategy.GREEDY,
    run_experiment=runner,
)
for r in history:
    ts = [d["t"] for d in r.design.designs]
    print(f"round {r.round}: k_hat = {r.estimation.parameters['k']:.3f}, batch t = {ts}")

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

round 0: k_hat = 0.220, batch t = [4.537760617533834, 4.537745350746717]
round 1: k_hat = 0.542, batch t = [1.844507759369574, 1.8445129653828807]
round 2: k_hat = 0.589, batch t = [1.6976164539462149, 1.697771094018572]

Notes and scope#

What is supported: joint D/A/E/ME optimality over \(N\) independent experiments with continuous design bounds, with an additive prior FIM carried forward by sequential_doe.

What is out of scope in v1:

  • Robust / expected-information designs \(\mathbb{E}_\theta[\text{FIM}]\) under parameter uncertainty (see Sandrin et al. [2025] for a modern framework).

  • Mixed-integer batch designs (choose \(N\) from a discrete candidate pool).

  • Bayesian-style q-acquisitions on GP surrogates; for that flavour of parallel experiment design see González and Zavala [2023].

  • Heterogeneous batches (different Experiment types in one batch).

  • End-to-end JAX autograd through the joint objective; the joint solver currently uses scipy finite differences.