Simulation#

Fake-data generators and process simulators: a bioreactor batch generator, a generator of data with a known latent structure, and a hidden-response-surface simulator an agent can probe one experiment at a time without being able to read the model behind it.

Bioreactor batch simulator#

A deterministic fed-batch bioreactor simulator with tunable disturbance channels.

Purpose#

This module is the quantitative baseline for batch trajectory adaptation and, later, batch mid-course correction. It makes three claims demonstrable, with numbers a reader can reproduce exactly from a seed:

Replaying a “golden batch” trajectory open-loop does not reproduce the golden outcome, because each batch carries its own disturbances.
The spread persists even when the measured initial conditions are held identical, because disturbances also arise during the batch.
The quality variance under replay therefore decomposes into a bucket that pre-batch (feedforward) adaptation can address, a bucket that only mid-course (feedback) correction can address, and a noise floor.

The model#

A 10-day fed-batch bioreactor with biomass X, substrate S, product (titer) P and working volume V. The specific growth rate follows the gamma-concept model of Rosso et al. (1995): a cardinal temperature model with inflection (CTMI) multiplied by a cardinal pH model (CPM), each equal to 1 at its optimum and exactly 0 outside its cardinal range, further multiplied by Monod substrate limitation, a within-batch disturbance phi(t), an initial-condition growth-inhibition factor and an oxygen-limitation factor:

mu_pot = mu_opt * gamma_T(T) * gamma_pH(pH) * S / (K_S + S) * phi * inh
f_O2   = smooth_min(1, our_capacity / ((o2_yield * mu_pot + o2_m) * X))
mu     = mu_pot * f_O2

k_d = k_d0 * (1 + exp((T - temp_death) / width_death))
      + k_d_starv * k_s_starv / (k_s_starv + S)
      + k_d_hyp * (1 - f_O2)
q_P = (alpha_lp * mu + beta_lp * gamma_q(T) * gamma_pH(pH) * phi * f_O2)
      * S / (k_sp + S)

dX/dt = mu * X - k_d * X - (F/V) * X
dS/dt = -(mu / yield_xs + maintenance) * X - q_P * X / yield_ps
        + (F/V) * (feed_substrate - S)
dP/dt = q_P * X - (F/V) * P
dV/dt = F

The couplings are chosen so the industrially standard biphasic temperature shift (a warm growth phase, then a mild-hypothermia production phase) is a genuine optimum rather than decoration, and so the right schedule depends on the batch’s own conditions:

The oxygen-transfer capacity caps the biomass pile the reactor can sustain, and the cap is temperature-dependent because warmer cells each demand more oxygen. Overshooting the cap is punished by hypoxic death, an irreversible loss, so the optimal pile size is interior.
Strong hypothermic growth arrest below about 28 degC (the reason biphasic mammalian-cell culture works) freezes the pile once the batch shifts cold, so the timing of the shift decides the production capacity, and the right timing depends on the inoculum and the growth rate the raw-material lot supports.
Production consumes substrate (yield_ps) and stalls when it runs out (k_sp); starvation kills cells (k_d_starv). Residual growth at too-warm production temperatures burns the feed that production needs, which is what places the production hold at an interior optimum (temp_production) below the isolated productivity optimum temp_q_opt.
The growth-inhibition factor inh (from the raw-material impurity latent factor) acts on growth only, so an inhibited lot needs a longer warm growth phase, not a scaled-down copy of the same batch.

The gamma hypothesis (temperature and pH acting independently on growth) is a modelling choice, not an established fact: published work finds the cardinal temperatures themselves shift with pH. It is adopted here because it keeps the true optimum interpretable in closed form.

Sensitivity by construction#

A simulator whose quality output moves visibly when an input moves by less than an instrument can resolve is not believable. Three structural properties keep this model insensitive to meaningless input changes while still responsive to real ones:

Quality is a time integral of bounded, smooth rates over 240 integration steps, so high-frequency input noise averages out and only sustained deviations accumulate.
The nominal recipe sits at stationary points of the response: pH is held at its cardinal optimum, and the production-phase temperature hold is the interior optimum of the hold-temperature response. Instrument-scale deviations around a stationary point are second order.
Sustained multi-degree deviations, by contrast, cross real mechanisms (hypothermic growth arrest on the cold side, feed-burning residual growth and hypoxia on the warm side), so they cost percent-level titer, as they should.

BioreactorSimulator.sensitivity_budget() computes the resulting budget from the live configuration, so the realism claim can be checked rather than taken on trust. Indicative values for the default configuration: zero-mean control-loop noise at instrument scale (sd 0.15 degC, 0.02 pH) moves the final titer by under 0.3% (standard deviation), a sustained 0.1 degC bias costs about 0.2%, a sustained 0.02 pH bias is invisible (about 0.01%), while a sustained 1 degC bias costs 7 to 11% and a sustained 2 degC bias about 20%. Around the operating point overheating costs more than the same excursion undercooling (feed burn plus hypoxia); far from it undercooling costs more (full growth arrest).

Disturbance channels#

Three independently tunable channels, each with a scale that can be set to zero (ic_scale, within_batch_scale, noise_scale):

Measured initial conditions: an 11-variable upstream Z block generated from three latent factors (seed viability, medium richness, inhibitor level) with three cluster centres, the feed classes A, B and C. Only the three latent factors drive the process, so Z also carries directions that do not matter, as real upstream data do.
Within-batch, unmeasured but observable: an Ornstein-Uhlenbeck process on log phi(t) with a correlation time comparable to the batch length, plus a per-batch feed-rate drift. It is not predictable from Z, is continuously visible in the oxygen and offgas trajectories, and its present partly predicts its own future, which is exactly why observing the running batch helps.
Control-loop and measurement noise: autocorrelated zero-mean deviation of the realised trajectory from its setpoints, plus independent measurement noise on every recorded tag, at instrument scale.

All random draws are made on every call and multiplied by their channel scale exactly once, so setting a scale to zero removes that channel without changing the draw sequence of the others: the same seed with a channel switched off is a true counterfactual for the same batch. The initial-condition scale acts where upstream data are drawn (sample_initial_conditions()); a caller-supplied Z block is measured data and is used at face value.

References

Rosso, L., Lobry, J. R., Bajard, S. and Flandrois, J. P. (1995), “Convenient model to describe the combined effects of temperature and pH on microbial growth”, Applied and Environmental Microbiology, 61 (2), 610-616.

Luedeking, R. and Piret, E. L. (1959), “A kinetic study of the lactic acid fermentation”, Journal of Biochemical and Microbiological Technology and Engineering, 1 (4), 393-412.

process_improve.simulation.batch.cardinal_temperature(temperature, t_min, t_opt, t_max)[source]#

Cardinal temperature model with inflection (CTMI) of Rosso et al. (1995).

Equal to 1 at t_opt, exactly 0 at and beyond the cardinal bounds, and smoothly curved in between with an inflection in the suboptimal range.

Parameters:

temperature (float or np.ndarray) – Temperature(s) [degC].
t_min (float) – Cardinal temperatures [degC]: no growth at or below t_min, maximum growth at t_opt, no growth at or above t_max. The CTMI is well posed only for t_opt >= (t_min + t_max) / 2; below that the denominator has a root strictly inside the cardinal window and the formula is meaningless there, so such triples are rejected.
t_opt (float) – Cardinal temperatures [degC]: no growth at or below t_min, maximum growth at t_opt, no growth at or above t_max. The CTMI is well posed only for t_opt >= (t_min + t_max) / 2; below that the denominator has a root strictly inside the cardinal window and the formula is meaningless there, so such triples are rejected.
t_max (float) – Cardinal temperatures [degC]: no growth at or below t_min, maximum growth at t_opt, no growth at or above t_max. The CTMI is well posed only for t_opt >= (t_min + t_max) / 2; below that the denominator has a root strictly inside the cardinal window and the formula is meaningless there, so such triples are rejected.

Returns:

The dimensionless growth factor gamma_T in [0, 1], with the same shape as temperature.

Return type:

np.ndarray

process_improve.simulation.batch.cardinal_ph(ph, ph_min, ph_opt, ph_max)[source]#

Cardinal pH model (CPM) of Rosso et al. (1995).

Equal to 1 at ph_opt and exactly 0 at and beyond the cardinal bounds.

Parameters:

ph (float or np.ndarray) – pH value(s) [-].
ph_min (float) – Cardinal pH values: no growth at or below ph_min, maximum growth at ph_opt, no growth at or above ph_max.
ph_opt (float) – Cardinal pH values: no growth at or below ph_min, maximum growth at ph_opt, no growth at or above ph_max.
ph_max (float) – Cardinal pH values: no growth at or below ph_min, maximum growth at ph_opt, no growth at or above ph_max.

Returns:

The dimensionless growth factor gamma_pH in [0, 1], with the same shape as ph.

Return type:

np.ndarray

class process_improve.simulation.batch.BioreactorConfig(mu_opt=0.8, k_s=0.2, yield_xs=0.4, maintenance=0.02, alpha_lp=0.3, beta_lp=0.35, yield_ps=0.9, k_sp=0.1, k_d_starv=0.08, k_s_starv=0.12, temp_min=27.5, temp_opt=36.8, temp_max=41.5, k_d0=0.015, temp_death=39.0, width_death=0.9, temp_q_min=22.0, temp_q_opt=31.5, temp_q_max=40.5, ph_min=6.3, ph_opt=7.1, ph_max=7.9, batch_days=10.0, samples_per_batch=20, steps_per_day=24, feed_rate=0.055, feed_substrate=50.0, volume_initial=1.0, biomass_initial=0.3, substrate_initial=5.0, temp_bounds=(28.0, 39.0), ph_bounds=(6.6, 7.6), shift_start_day=3.0, shift_end_day=4.5, temp_production=29.05, ic_scale=1.0, within_batch_scale=1.0, noise_scale=1.0, ou_tau_days=6.0, ou_sd=0.15, feed_drift_sd=0.06, control_sd_temp=0.15, control_sd_ph=0.02, control_tau_h=6.0, meas_sd_temp=0.05, meas_sd_ph=0.01, meas_sd_do=1.0, meas_sd_co2=0.08, meas_sd_volume=0.005, o2_yield=0.55, o2_maintenance=0.02, our_capacity=1.85, k_d_hyp=0.4, rq=1.0, co2_gain=1.0, co2_inlet_pct=0.04)[source]#

Bases: object

Kinetic, operating and disturbance parameters of the simulated bioreactor.

All parameters have defaults chosen so the nominal batch lands in a physiologically plausible range for a pilot-scale mammalian-cell fed-batch process (peak biomass a few g/L, final titer of order 10 g/L). Units are given per field. The dataclass is frozen; derive variants with dataclasses.replace().

Parameters:

mu_opt (float) – Maximum specific growth rate at the cardinal optima [1/day].
k_s (float) – Monod half-saturation constant for the substrate [g/L].
yield_xs (float) – Biomass yield on substrate [g biomass / g substrate].
maintenance (float) – Maintenance substrate consumption [g substrate / g biomass / day].
alpha_lp (float) – Luedeking-Piret coefficients: growth-associated product yield [g product / g biomass] and non-growth-associated specific productivity [g product / g biomass / day].
beta_lp (float) – Luedeking-Piret coefficients: growth-associated product yield [g product / g biomass] and non-growth-associated specific productivity [g product / g biomass / day].
yield_ps (float) – Product yield on substrate [g product / g substrate]: production consumes substrate, so a large biomass pile competes with its own productivity for feed.
k_sp (float) – Half-saturation substrate concentration for product formation [g/L]: production stalls when the substrate runs out.
k_d_starv (float) – Additional death rate under full starvation [1/day].
k_s_starv (float) – Substrate concentration at which the starvation death rate reaches half its maximum [g/L].
temp_min (float) – Cardinal temperatures for growth [degC].
temp_opt (float) – Cardinal temperatures for growth [degC].
temp_max (float) – Cardinal temperatures for growth [degC].
k_d0 (float) – Baseline first-order death rate [1/day].
temp_death (float) – Onset temperature [degC] and width [degC] of the exponential rise in the death rate.
width_death (float) – Onset temperature [degC] and width [degC] of the exponential rise in the death rate.
temp_q_min (float) – Cardinal temperatures for the non-growth-associated productivity [degC]; temp_q_opt sits below temp_opt, which is what makes the biphasic temperature shift a real trade-off.
temp_q_opt (float) – Cardinal temperatures for the non-growth-associated productivity [degC]; temp_q_opt sits below temp_opt, which is what makes the biphasic temperature shift a real trade-off.
temp_q_max (float) – Cardinal temperatures for the non-growth-associated productivity [degC]; temp_q_opt sits below temp_opt, which is what makes the biphasic temperature shift a real trade-off.
ph_min (float) – Cardinal pH values, shared by growth and productivity.
ph_opt (float) – Cardinal pH values, shared by growth and productivity.
ph_max (float) – Cardinal pH values, shared by growth and productivity.
batch_days (float) – Batch duration [day].
samples_per_batch (int) – Number of recorded samples (and setpoint intervals) per batch.
steps_per_day (int) – Integration steps per day for the fixed-step RK4 integrator.
feed_rate (float) – Nominal constant feed rate [L/day].
feed_substrate (float) – Substrate concentration in the feed [g/L].
volume_initial (float) – Working volume at inoculation [L].
biomass_initial (float) – Nominal initial biomass and substrate concentrations [g/L]; the measured initial-condition channel perturbs these per batch.
substrate_initial (float) – Nominal initial biomass and substrate concentrations [g/L]; the measured initial-condition channel perturbs these per batch.
temp_bounds (tuple[float, float]) – Recipe-allowed operating window for the setpoints; requested and realised trajectories are validated or clipped against these.
ph_bounds (tuple[float, float]) – Recipe-allowed operating window for the setpoints; requested and realised trajectories are validated or clipped against these.
shift_start_day (float) – Start and end [day] of the nominal biphasic temperature ramp from temp_opt down to temp_production.
shift_end_day (float) – Start and end [day] of the nominal biphasic temperature ramp from temp_opt down to temp_production.
temp_production (float) – The production-phase hold temperature of the nominal recipe [degC]. It sits below the isolated productivity optimum temp_q_opt because residual growth at warmer temperatures consumes the feed that production needs; cold enough to arrest growth is what a sensible recipe holds.
ic_scale (float) – Scale of the upstream variation put into drawn initial conditions (campaigns that draw their own Z block, and sample_initial_conditions() when called through them); 0 collapses drawn batches to the nominal upstream values. A caller-supplied Z block is measured data and is always used at face value, independent of this scale.
within_batch_scale (float) – Scale of the unmeasured within-batch channel; 0 switches it off.
noise_scale (float) – Scale of the control-loop and measurement noise channel; 0 switches it off.
ou_tau_days (float) – Correlation time of the Ornstein-Uhlenbeck disturbance on log phi(t) [day].
ou_sd (float) – Stationary standard deviation of log phi(t) [-].
feed_drift_sd (float) – Per-batch fractional standard deviation of the realised feed rate [-].
control_sd_temp (float) – Stationary standard deviation of the realised-minus-setpoint control error for temperature [degC] and pH [-].
control_sd_ph (float) – Stationary standard deviation of the realised-minus-setpoint control error for temperature [degC] and pH [-].
control_tau_h (float) – Correlation time of the control error [hour].
meas_sd_temp (float) – Measurement noise standard deviations on the recorded tags [degC, -, % saturation, % offgas, L].
meas_sd_ph (float) – Measurement noise standard deviations on the recorded tags [degC, -, % saturation, % offgas, L].
meas_sd_do (float) – Measurement noise standard deviations on the recorded tags [degC, -, % saturation, % offgas, L].
meas_sd_co2 (float) – Measurement noise standard deviations on the recorded tags [degC, -, % saturation, % offgas, L].
meas_sd_volume (float) – Measurement noise standard deviations on the recorded tags [degC, -, % saturation, % offgas, L].
o2_yield (float) – Oxygen demand per unit growth [g O2 / g biomass].
o2_maintenance (float) – Maintenance oxygen uptake [g O2 / g biomass / day].
our_capacity (float) – Maximum oxygen transfer the sparger and agitation can supply [g O2 / L / day]. This is the binding constraint of the process: it caps the biomass pile the reactor can sustain, and the cap is temperature-dependent because warmer cells demand more oxygen each.
k_d_hyp (float) – Additional death rate under full oxygen limitation [1/day]: overshooting the sustainable biomass pile is an irreversible loss, which is what gives the optimal schedule curvature on both sides.
rq (float) – Respiratory quotient, the CO2 evolved per O2 consumed [g/g].
co2_gain (float) – Offgas CO2 percentage per unit CO2 evolution rate [% / (g/L/day)].
co2_inlet_pct (float) – CO2 percentage of the inlet gas [%].

property n_steps: int#: Total number of RK4 integration steps per batch.

property interval_days: float#: Duration of one setpoint interval (one recorded sample) [day].

property sample_days: ndarray#: Recording times, the end of each setpoint interval [day].

property interval_start_days: ndarray#: Setpoint interval start times, the index of a trajectory frame [day].

process_improve.simulation.batch.sample_initial_conditions(n_batches, *, proportions=None, ic_scale=1.0, random_state=None)[source]#

Draw upstream (initial-condition) data for a campaign of batches.

Each batch belongs to one of three feed-disturbance classes, A, B or C, with a cluster centre in the three-dimensional latent space (seed viability, medium richness, inhibitor level). The observed Z block has 11 variables generated from those latent factors through a fixed loading matrix plus residual noise, so a PCA of Z recovers about three meaningful components and shows the three classes as clusters.

Parameters:

n_batches (int) – Number of batches to draw.
proportions (dict[str, float], optional) – Expected class proportions keyed by "A", "B", "C". They must be non-negative and sum to a positive number (they are normalised). Default {"A": 0.40, "B": 0.35, "C": 0.25}.
ic_scale (float, default=1.0) – Scale of the initial-condition variation. 0 collapses every batch to the nominal upstream values while keeping the class labels.
random_state (int, np.random.Generator, or None) – Seed or generator; see process_improve._random.check_random_state.

Returns:

result – With keys z (DataFrame, n_batches rows by 11 upstream variables, integer batch ids as index), classes (Series of “A” / “B” / “C” labels, same index) and latent (DataFrame of the three latent factor values, same index).

Return type:

sklearn.utils.Bunch

class process_improve.simulation.batch.BioreactorSimulator(config=None)[source]#

Bases: object

A deterministic fed-batch bioreactor with tunable disturbance channels.

See the module docstring for the model, its citations, and the design constraints. Given the same configuration, inputs and random_state, every method reproduces its results exactly.

Parameters:: config (BioreactorConfig, optional) – The full parameter set; defaults to BioreactorConfig().

Examples

>>> from process_improve.simulation import BioreactorSimulator
>>> sim = BioreactorSimulator()
>>> golden = sim.golden_trajectory()
>>> campaign = sim.simulate_campaign(50, policy="replay", trajectory=golden.trajectory, random_state=42)
>>> float(campaign.quality["titer"].std()) > 0.0
True

nominal_trajectory()[source]#

Return the nominal (recipe) setpoint schedule: biphasic temperature, pH held.

Temperature holds the growth optimum until shift_start_day, ramps linearly to the production hold temp_production by shift_end_day, and holds it for the rest of the batch. pH is held at its cardinal optimum throughout, matching the industrial practice of holding pH and shifting temperature.

Returns:: samples_per_batch rows indexed by setpoint interval start time [day], columns ["pH", "temperature"]. Each row is the setpoint held over the following interval (zero-order hold).
Return type:: pd.DataFrame

simulate_batch(initial_conditions=None, trajectory=None, *, random_state=None)[source]#

Simulate one batch under given initial conditions and setpoints.

Parameters:

initial_conditions (pd.Series, optional) – One row of the upstream Z block (index: the 11 upstream variable names). None runs the nominal batch (all upstream variables at their means).
trajectory (pd.DataFrame, optional) – Setpoint schedule with samples_per_batch rows and columns ["pH", "temperature"]; each row is held over its interval (zero-order hold). None uses nominal_trajectory().
random_state (int, np.random.Generator, or None) – Seed or generator for the disturbance and noise draws.

Returns:

result – With keys:

tags: DataFrame, samples_per_batch rows indexed by sample time [day], columns ["pH", "temperature", "dissolved_oxygen", "offgas_co2", "volume"]: what the historian records, including measurement noise. Each row reads at the end of its interval: states at that instant, inputs as delivered over the interval that ends there.
titer: float, the final product concentration [g/L].
states: DataFrame on the integration grid (indexed by day) with columns ["biomass", "substrate", "titer", "volume", "phi", "temperature", "pH"]: the noise-free god view, including the unmeasured disturbance phi. The input columns (phi, temperature, pH) hold the value in force during the step starting at each grid instant.
realised_trajectory: DataFrame like tags but only ["pH", "temperature"] and without measurement noise: what the control loops delivered over the interval ending at each sample time.
initial_conditions: Series, the Z row used, taken at face value (see ic_scale in BioreactorConfig).

Return type:

sklearn.utils.Bunch

optimal_trajectory(initial_conditions=None, *, n_knots=4, n_starts=5, random_state=0)[source]#

Optimise the setpoint schedule for a batch’s initial conditions on the simulator’s own model.

Maximises the deterministic (disturbance-free) final titer over pH and temperature schedules parameterised by n_knots values each, linearly interpolated across the batch and held per interval. Because the optimiser queries the simulator’s own model rather than a fitted one, the result is the feedforward reference any data-driven scheme can be scored against. The search is a multi-start local one within the knot parameterisation, so its titer is a lower bound on the true optimum for these initial conditions.

Parameters:

initial_conditions (pd.Series, optional) – One upstream Z row; None for the nominal batch.
n_knots (int, default=4) – Number of knots per manipulated variable.
n_starts (int, default=5) – Multi-start count: the nominal trajectory plus n_starts - 1 random starts inside the operating bounds.
random_state (int, np.random.Generator, or None, default=0) – Seed for the random starts; the default makes the search reproducible without an argument.

Returns:

result – With keys trajectory (DataFrame in the same layout as nominal_trajectory()), titer (float, the deterministic titer of that trajectory) and optimizer_success (bool).

Return type:

sklearn.utils.Bunch

golden_trajectory(*, n_knots=4, n_starts=5, random_state=0)[source]#

Find the golden batch: the optimal schedule for the nominal initial conditions.

This is what golden-batch practice enshrines as the recipe. It is the true optimum only for the conditions under which it was found; the whole point of the baseline is what happens when it is replayed under other conditions. Parameters and return value match optimal_trajectory() with initial_conditions=None.

Parameters:

n_knots (int)
n_starts (int)
random_state (int | Generator | None)

Return type:

Bunch

simulate_campaign(n_batches, *, policy='replay', trajectory=None, initial_conditions=None, mv_variation=0.0, n_knots=4, n_starts=5, random_state=None)[source]#

Simulate a campaign of batches under a named operating policy.

Parameters:

n_batches (int) – Number of batches in the campaign.
policy ({"replay", "historical", "adapted"}, default="replay") –
- "replay": every batch is given the same setpoint schedule, the golden-batch practice.
- "historical": the same schedule plus deliberate per-batch setpoint variation of size mv_variation (see below). A perfectly consistent history contains no information about how the controls affect quality; this policy produces a history that does.
- "adapted": each batch runs the schedule that maximises the disturbance-free titer for its own initial conditions, computed from the simulator’s model via optimal_trajectory(). This estimates the ceiling a feedforward scheme with a perfect process model could reach; it is not an implementable policy.
trajectory (pd.DataFrame, optional) – The schedule used by "replay" and "historical"; defaults to nominal_trajectory(). Ignored by "adapted".
initial_conditions (pd.DataFrame, optional) – The upstream Z block, one row per batch (the 11 upstream variables as columns). None draws it with sample_initial_conditions() using the configuration’s ic_scale.
mv_variation (float, default=0.0) – Size of the deliberate variation for "historical": each batch draws independent offsets at five evenly spaced knots, linearly interpolated over the batch, for temperature (standard deviation mv_variation degC per knot) and for pH (standard deviation 0.1 * mv_variation per knot), clipped to the operating bounds. The knot basis matters for identification: it excites local schedule shapes (an early-only or late-only move), which a constant-plus-ramp variation cannot, so a regression fitted on the campaign carries causal information about partial-batch moves. Flores-Cerrillo and MacGregor (2004) identified their mid-course-correction model from batches with deliberate moves in the manipulated trajectories at the decision points, and suggest batch-to-batch updating for a model identified from historical data alone.
n_knots (int) – Passed to optimal_trajectory() for the "adapted" policy, which runs the optimiser once per batch; lower values trade optimality for speed. Ignored by the other policies.
n_starts (int) – Passed to optimal_trajectory() for the "adapted" policy, which runs the optimiser once per batch; lower values trade optimality for speed. Ignored by the other policies.
random_state (int, np.random.Generator, or None) – Seed or generator; child seeds are spawned per batch, so a campaign is reproducible end to end.

Returns:

result – With keys:

batches: dict[int, pd.DataFrame], the recorded tags per batch in the package’s standard batch-dictionary format.
quality: DataFrame indexed by batch id with the single column titer [g/L].
initial_conditions: DataFrame, the Z block used.
classes: Series of feed-class labels (“A”/”B”/”C”), or "?" when initial_conditions was supplied by the caller.
trajectories: dict[int, pd.DataFrame], the requested setpoint schedule per batch.

Return type:

sklearn.utils.Bunch

sensitivity_budget(*, n_noise_replicates=100, random_state=0)[source]#

How much the final titer moves under standard input perturbations.

This is the realism check: a credible process model must not respond visibly to input changes smaller than an instrument can resolve, and must respond clearly to sustained multi-degree deviations. The rows cover zero-mean control-loop noise at instrument scale, sustained setpoint biases of increasing size, and a single-sample excursion. All effects are measured on the deterministic model (disturbance channels off) except the control-noise row, which uses the configured control-error model.

Parameters:

n_noise_replicates (int, default=100) – Number of batches used for the control-loop-noise row.
random_state (int, np.random.Generator, or None, default=0) – Seed for the control-noise replicates.

Returns:

One row per perturbation with columns perturbation (index), titer_g_L and effect_pct (percentage change of the final titer against the unperturbed nominal batch; for the noise row, the standard deviation across replicates).

Return type:

pd.DataFrame

process_improve.simulation.batch.variance_decomposition(simulator, n_batches=200, *, trajectory=None, random_state=None)[source]#

Split the replay-policy titer variance into its three sources.

Runs four replay campaigns with the same size and schedule: all channels on (the total), then each channel alone. The initial-condition bucket is what pre-batch (feedforward) trajectory adaptation could remove; the within-batch bucket is what only mid-course (feedback) correction can reach; the noise bucket is the floor. Because the model is nonlinear the three buckets do not sum exactly to the total; the difference is reported as the interaction residual rather than hidden.

Parameters:

simulator (BioreactorSimulator) – The simulator whose configuration (including channel scales) defines the “all channels on” case.
n_batches (int, default=200) – Batches per campaign.
trajectory (pd.DataFrame, optional) – The replayed schedule; defaults to the simulator’s nominal.
random_state (int, np.random.Generator, or None) – Seed or generator; each campaign draws its own child seed.

Returns:

Rows ["measured initial conditions", "within-batch disturbance", "control and measurement noise", "interaction residual", "total"] with columns variance [g^2/L^2], sd [g/L], cv_pct (sd as a percentage of the all-channels mean titer) and pct_of_total (share of the total variance; the residual’s share can be negative).

Return type:

pd.DataFrame

Known latent structure#

Kevin Dunn, 2010-2026. MIT License.

Data with a known latent structure, for checking latent-variable methods.

LatentStructure describes a linear process

\[\mathbf{X} = \mathbf{T}\mathbf{P}^\top + \mathbf{E}, \qquad \mathbf{Y} = \mathbf{T}\mathbf{B} + \mathbf{F},\]

in which the caller sets the number of latent variables, how strongly each one varies, and which of them drive each response. LatentStructure.sample() draws data from it, optionally with cells missing completely at random. A second call with another seed draws fresh rows from the same process, which is the ground truth a cross-validated estimate can be checked against. The best that any linear model can do on new rows is known in closed form (LatentStructure.population_coefficients() and LatentStructure.population_r2()).

The loadings are columns of a Hadamard matrix scaled to unit length, so every entry is \(\pm 1/\sqrt{K}\). Every X column then has the same variance, and autoscaling, which PCA and PLS apply by default, is a uniform rescale that leaves the structure unchanged. With loadings of unequal row norms, autoscaling makes the noise heteroscedastic, and a PLS model of the autoscaled data then needs more components than there are latent variables to reach the best linear predictor, so the number of latent variables the process was built with would no longer be the right answer.

class process_improve.simulation.latent.LatentStructure(x_sd, y_coefficients, n_features=16, noise_sd=0.3, y_noise_sd=None)[source]#

Bases: object

A linear latent-variable process whose structure is known.

Parameters:

x_sd (sequence of float) – Standard deviation of each latent variable. Its length is the number of latent variables in X.
y_coefficients (array-like of shape (n_latent,) or (n_latent, n_targets)) – B: how each latent variable drives each response. A zero row is variation in X that Y does not see (Y-orthogonal variation). A 1-D array means one response.
n_features (int, default=16) – Number of X variables, K: a power of two, larger than the number of latent variables.
noise_sd (float, default=0.3) – Standard deviation of the noise added to every X cell.
y_noise_sd (float, optional) – Standard deviation of the noise added to every Y cell; noise_sd when omitted.

Examples

Two latent variables drive y; a PLS model should need two components.

>>> process = LatentStructure(x_sd=[3.0, 1.0], y_coefficients=[1.0, 1.0])
>>> train = process.sample(60, random_state=0)
>>> fresh = process.sample(4000, random_state=1)  # new rows from the same process
>>> process.n_relevant
2

property n_latent: int#: Number of latent variables in X.

property n_targets: int#: Number of responses in Y.

property n_relevant: int#: Number of latent variables that drive at least one response (a nonzero row of B).

property loadings: DataFrame#

orthonormal columns with every entry equal to +-1/sqrt(K).

Type:: The loadings P

population_coefficients()[source]#

Coefficients of the best linear predictor of Y from X, on an infinite sample.

With \(\boldsymbol{\Sigma}_{XX} = \mathbf{P}\boldsymbol{\Lambda}\mathbf{P}^\top + \sigma^2\mathbf{I}\) and orthonormal \(\mathbf{P}\), the regression \(\boldsymbol{\Sigma}_{XX}^{-1}\boldsymbol{\Sigma}_{XY}\) reduces to \(\mathbf{P}\,\mathrm{diag}\{\lambda_j / (\lambda_j + \sigma^2)\}\,\mathbf{B}\): the noise in X shrinks each latent variable’s contribution by the fraction of its variance that X reveals. Every variable has mean zero, so there is no intercept.

Return type:: pd.DataFrame of shape (n_features, n_targets)

population_r2()[source]#

Fraction of each response’s variance that the best linear predictor explains on new rows.

No linear model does better on new rows from this process, so this is the value a cross-validated \(Q^2\) estimates at best (an estimate from a finite sample can land slightly above it). The predictable part is \(\sum_j b_{jm}^2\,\lambda_j^2 / (\lambda_j + \sigma^2)\), out of a total variance \(\sum_j b_{jm}^2\,\lambda_j + \sigma_y^2\).

Return type:: pd.Series of length n_targets

sample(n_samples, *, missing_fraction=0.0, random_state=None)[source]#

Draw n_samples rows from the process.

The scores, the X noise and the Y noise are drawn in that order, and the missing cells after them, so the same random_state gives the same values with or without missing cells.

Parameters:

n_samples (int) – Number of rows, N.
missing_fraction (float, default=0.0) – Fraction of the X and of the Y cells deleted completely at random (set to NaN), in [0, 1). A row never loses every X cell: one is restored at random, so every row can still be scored.
random_state (int, numpy.random.Generator or None, default=None) – Seed or generator; see process_improve._random.check_random_state.

Returns:

With X (DataFrame, N x K), Y (DataFrame, N x M) and scores (DataFrame, N x n_latent, the true latent values T).

Return type:

sklearn.utils.Bunch

Hidden response-surface simulator#

Kevin Dunn, 2010-2026. MIT License.

Deterministic fake-data simulator for DOE demonstrations.

The simulator is fully specified by a JSON-serialisable private_state dict containing:

seed: int - base RNG seed for coefficient generation.
factors: list of {name, low, high, units?}.
outputs: list of {name, units?, direction?}.
structural_hints: list of free-text hints (“negative interaction between pH and surfactant”, etc).
noise_level: "low" | "medium" | "high".
time_drift: bool.
model_version: int (schema version for future migrations).

Given this state, materialize_model() regenerates the exact same response-surface coefficients on every call, across processes, across machines. simulate() evaluates that surface at a given factor setting and adds fresh Gaussian noise each call so identical inputs yield similar-but-not-identical outputs - matching the behaviour of a real physical asset.

process_improve.simulation.model.validate_factors(factors)[source]#

Validate the factors list from a create_simulator call.

Parameters:: factors (list[dict[str, Any]])
Return type:: None

process_improve.simulation.model.validate_outputs(outputs)[source]#

Validate the outputs list from a create_simulator call.

Parameters:: outputs (list[dict[str, Any]])
Return type:: None

process_improve.simulation.model.validate_noise_level(noise_level)[source]#

Reject any noise_level outside the fixed enum.

Parameters:: noise_level (str)
Return type:: None

process_improve.simulation.model.materialize_model(private_state)[source]#

Regenerate the full coefficient set from a private_state dict.

Deterministic: identical private_state always returns identical coefficients, across processes and machines. This is the secret the simulator hides from the LLM - callers that have private_state already have the model, so “revealing” is free.

Parameters:: private_state (dict[str, Any])
Return type:: dict[str, Any]

process_improve.simulation.model.simulate(private_state, settings, timestamp_offset_days=0.0, *, random_state=None)[source]#

Evaluate the hidden response surface at settings, with fresh noise.

Parameters:

private_state (dict[str, Any]) – The state dict produced by create_simulator (and persisted by the host). Passed in whole so this function stays pure.
settings (dict[str, float]) – Mapping of factor-name → numeric value, in the same units as the factor’s declared range. Missing factors are filled with the mid-range value (with a warning). Out-of-range values are clipped to the declared bounds (with a warning).
timestamp_offset_days (float) – Optional time axis passed by the caller when time_drift is enabled on the simulator. Ignored otherwise.
random_state (int | Generator | None) – Seeds the measurement noise, per the reproducibility contract. The default, None, draws fresh entropy on every call, which is what a simulator standing in for a real process should do: two runs at the same settings must not return the same number. Pass an int or a numpy.random.Generator when a run has to be repeatable, for example to pin a test or to hand a reader a worked example they can reproduce. It is deliberately not part of the simulate_process tool contract, so a model driving the simulator cannot freeze its noise.

Returns:

{settings, outputs, warnings, timestamp_offset_days}.

Return type:

dict

process_improve.simulation.model.draw_initial_seed()[source]#

Return a non-negative int suitable for use as a seed in private_state.

Uses secrets.randbits(63) to draw 63 bits of cryptographic entropy. 63 (not 64) keeps the value in a positive signed 64-bit int so JSON parsers and SQLAlchemy Integer columns can round-trip it without width worries. Previous versions truncated SeedSequence().entropy to 31 bits, which gave an attacker who could observe simulator outputs a brute-forceable seed space (SEC-28 / #277). With 63 bits, enumerating the seed space is no longer feasible.

Return type:: int

Host-injected state#

Kevin Dunn, 2010-2026. MIT License.

Host-side injection channel for the fake-data simulator tools.

simulate_process and reveal_simulator need two pieces of state that the LLM must never be able to supply itself:

simulator_state - the hidden response-surface model the host persists per conversation.
confirmed - whether the user has cleared the double-confirmation gate guarding a model reveal.

Earlier these travelled as keyword arguments deliberately kept out of the public JSON schema, on the assumption the host would inject them server-side. But the dispatch path forwards the tool input verbatim as **kwargs, so a prompt-injected agent could re-introduce them as ordinary kwargs and bypass the reveal gate (SEC-15). The schema filter in process_improve.tool_spec.execute_tool_call() now drops such keys, and this module removes them from the kwarg surface entirely: they live in contextvars.ContextVar slots that only the host populates, so they cannot be reached through a tool call even if the filter regresses.

Host usage#

The host wraps each dispatch in simulator_host_context():

from process_improve.simulation.context import simulator_host_context
from process_improve.tool_spec import execute_tool_call

with simulator_host_context(simulator_state=state, confirmed=user_ok):
    result = execute_tool_call(tool_name, tool_input)

Outside such a block the context is empty, so simulate_process reports simulator_state_missing and reveal_simulator returns confirmation_needed - the safe defaults.

process_improve.simulation.context.get_injected_simulator_state()[source]#

Return the host-injected simulator_state for the current dispatch.

Returns None when called outside a simulator_host_context() block, which the simulator tools surface as simulator_state_missing.

Return type:: dict[str, Any] | None

process_improve.simulation.context.get_reveal_confirmed()[source]#

Return whether the host has confirmed a reveal for the current dispatch.

Defaults to False outside a simulator_host_context() block so an un-gated call can never expose the hidden model.

Return type:: bool

process_improve.simulation.context.simulator_host_context(*, simulator_state=None, confirmed=False)[source]#

Inject simulator state / reveal confirmation for the enclosed dispatch.

Parameters:

simulator_state (dict[str, Any] | None) – The hidden model the host persists for this conversation, or None to leave it unavailable.
confirmed (bool) – True only after the user has cleared the reveal double-confirmation gate. Coerced to bool.

Return type:

Iterator[None]

Notes

The values are stored in contextvars.ContextVar slots and reset on exit, so concurrent hosts (e.g. async tasks) do not see each other’s state. ContextVar values do not cross process boundaries, so a host that dispatches through the subprocess-isolated process_improve.tool_safety.safe_execute_tool_call() will not leak state into the worker; that untrusted-transport path is expected to report simulator_state_missing rather than carry persisted state.

Recipes and agent-callable tools#

Kevin Dunn, 2010-2026. MIT License.

Analysis recipes for the batch (bioreactor) simulator.

One guided workflow: establish, quantitatively, why replaying a golden batch’s schedule does not reproduce its outcome, and how the remaining variance splits between what can be addressed before a batch and what can only be addressed while it runs. Registered into the package-wide catalog on import; see process_improve.recipes for the framework.

Kevin Dunn, 2010-2026. MIT License.

Agent-callable tool wrappers for the simulation subpackage.

The tools exposed here:

create_simulator - records a hidden response-surface model.
simulate_process - evaluates the hidden model at given factor settings, with fresh Gaussian noise each call.
reveal_simulator - returns the underlying coefficients, gated behind a confirmed flag enforced by the host.
simulate_batch_campaign - runs a fed-batch bioreactor campaign under a named operating policy and reports the titer outcomes.
decompose_batch_quality_variance - splits the replay-campaign titer variance into its disturbance sources.
correct_batch_midcourse - corrects a handful of fresh batches at a decision point and reports the executed (re-simulated) outcomes.
evaluate_batch_control_policy - the campaign-level executed comparison of replay vs mid-course correction (optionally with the feedforward and oracle ceilings).

The three DOE tools above are stateless in the private_state sense described below; the four bioreactor tools have no hidden state at all (their model is deliberately open, see process_improve.simulation.batch) and are reproducible from their random_state field alone.

The tools are intentionally stateless: the hidden model lives in a private_state dict that the host (e.g. the factorial web app) persists and injects on each call. simulator_state and the reveal confirmed flag are not kwargs: the host supplies them out of band through process_improve.simulation.context.simulator_host_context(), which stores them in contextvars.ContextVar slots. Keeping them off the kwarg surface means a prompt-injected agent cannot re-introduce them through the dispatch path to forge state or bypass the reveal gate (SEC-15).

Pydantic input contract (ENG-04 / ENG-10): each tool pairs its @tool_spec decorator with a BaseModel carrying ConfigDict(extra="forbid"); the function receives the parsed model as its single positional argument. create_simulator.seed is declared as SkipJsonSchema[int | None] so the field is callable from Python (tests, notebooks) but does NOT appear in the JSON Schema exposed to the LLM.

class process_improve.simulation.tools.CreateSimulatorInput(*, process_description, factors, outputs, structural_hints=None, noise_level='medium', time_drift=False, seed=None)[source]#

Bases: BaseModel

Input contract for create_simulator.

Parameters:

process_description (Annotated[str, MinLen(min_length=1)])
factors (Annotated[list[dict[str, Any]], MinLen(min_length=1)])
outputs (Annotated[list[dict[str, Any]], MinLen(min_length=1)])
structural_hints (list[str] | None)
noise_level (Literal['low', 'medium', 'high'])
time_drift (bool)
seed (Annotated[int | None, SkipJsonSchema()])

model_config = {'extra': 'forbid'}#: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

process_improve.simulation.tools.create_simulator(spec)[source]#

Create a simulator; see tool spec for parameter details.

Parameters:: spec (CreateSimulatorInput)
Return type:: dict

class process_improve.simulation.tools.SimulateProcessInput(*, sim_id, settings, timestamp_offset_days=0.0)[source]#

Bases: BaseModel

Input contract for simulate_process.

Parameters:

sim_id (str)
settings (dict[str, float])
timestamp_offset_days (float)

model_config = {'extra': 'forbid'}#: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

process_improve.simulation.tools.simulate_process(spec)[source]#

Evaluate a simulator at settings.

simulator_state is not a parameter: the host injects it out of band via process_improve.simulation.context.simulator_host_context(), so it can never arrive as an LLM-supplied kwarg. If it is missing we return a structured error rather than dispatch with a guessed state.

Parameters:: spec (SimulateProcessInput)
Return type:: dict

class process_improve.simulation.tools.RevealSimulatorInput(*, sim_id)[source]#

Bases: BaseModel

Input contract for reveal_simulator.

Parameters:: sim_id (str)

model_config = {'extra': 'forbid'}#: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

process_improve.simulation.tools.reveal_simulator(spec)[source]#

Return the materialised model.

simulator_state and the confirmed flag are injected by the host out of band via process_improve.simulation.context.simulator_host_context(), never by the LLM. They are not parameters, so a tool call cannot forge state or pre-clear the double-confirmation gate.

Parameters:: spec (RevealSimulatorInput)
Return type:: dict

class process_improve.simulation.tools.SimulateBatchCampaignInput(*, n_batches, policy='replay', mv_variation=0.0, ic_scale=1.0, within_batch_scale=1.0, noise_scale=1.0, random_state=0)[source]#

Bases: BaseModel

Input contract for simulate_batch_campaign.

Parameters:

n_batches (Annotated[int, Ge(ge=2), Le(le=500)])
policy (Literal['replay', 'historical'])
mv_variation (Annotated[float, Ge(ge=0.0), Le(le=3.0)])
ic_scale (Annotated[float, Ge(ge=0.0), Le(le=3.0)])
within_batch_scale (Annotated[float, Ge(ge=0.0), Le(le=3.0)])
noise_scale (Annotated[float, Ge(ge=0.0), Le(le=3.0)])
random_state (Annotated[int, Ge(ge=0), Le(le=2147483647)])

model_config = {'extra': 'forbid'}#: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

process_improve.simulation.tools.simulate_batch_campaign(spec)[source]#

Simulate a bioreactor campaign and summarise the final-quality outcomes.

Parameters:: spec (SimulateBatchCampaignInput)
Return type:: dict[str, Any]

class process_improve.simulation.tools.DecomposeBatchQualityVarianceInput(*, n_batches=150, ic_scale=1.0, within_batch_scale=1.0, noise_scale=1.0, random_state=0)[source]#

Bases: BaseModel

Input contract for decompose_batch_quality_variance.

Parameters:

n_batches (Annotated[int, Ge(ge=10), Le(le=400)])
ic_scale (Annotated[float, Ge(ge=0.0), Le(le=3.0)])
within_batch_scale (Annotated[float, Ge(ge=0.0), Le(le=3.0)])
noise_scale (Annotated[float, Ge(ge=0.0), Le(le=3.0)])
random_state (Annotated[int, Ge(ge=0), Le(le=2147483647)])

model_config = {'extra': 'forbid'}#: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

process_improve.simulation.tools.decompose_batch_quality_variance(spec)[source]#

Decompose replay-campaign titer variance into its disturbance sources.

Parameters:: spec (DecomposeBatchQualityVarianceInput)
Return type:: dict[str, Any]

class process_improve.simulation.tools.CorrectBatchMidcourseInput(*, y_target=8.0, n_batches=5, decision_point=8, n_train=200, mv_variation=2.5, dead_band=1.0, random_state=0)[source]#

Bases: BaseModel

Input contract for the single-decision-point correction demonstration.

Parameters:

y_target (Annotated[float, Gt(gt=0.0), Le(le=50.0)])
n_batches (Annotated[int, Ge(ge=1), Le(le=10)])
decision_point (Annotated[int, Ge(ge=1), Le(le=19)])
n_train (Annotated[int, Ge(ge=30), Le(le=400)])
mv_variation (Annotated[float, Ge(ge=0.1), Le(le=5.0)])
dead_band (Annotated[float, Ge(ge=0.0), Le(le=10.0)])
random_state (Annotated[int, Ge(ge=0), Le(le=2147483647)])

model_config = {'extra': 'forbid'}#: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

process_improve.simulation.tools.correct_batch_midcourse(spec)[source]#

Correct fresh batches at one decision point and report the executed outcomes.

Parameters:: spec (CorrectBatchMidcourseInput)
Return type:: dict[str, Any]

class process_improve.simulation.tools.EvaluateBatchControlPolicyInput(*, y_target=8.0, n_train=200, n_test=40, mv_variation=2.5, decision_point=8, dead_band=1.0, include_ceilings=False, random_state=0)[source]#

Bases: BaseModel

Input contract for the campaign-level executed policy comparison.

Parameters:

y_target (Annotated[float, Gt(gt=0.0), Le(le=50.0)])
n_train (Annotated[int, Ge(ge=30), Le(le=400)])
n_test (Annotated[int, Ge(ge=5), Le(le=60)])
mv_variation (Annotated[float, Ge(ge=0.1), Le(le=5.0)])
decision_point (Annotated[int, Ge(ge=1), Le(le=19)])
dead_band (Annotated[float, Ge(ge=0.0), Le(le=10.0)])
include_ceilings (bool)
random_state (Annotated[int, Ge(ge=0), Le(le=2147483647)])

model_config = {'extra': 'forbid'}#: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

process_improve.simulation.tools.evaluate_batch_control_policy(spec)[source]#

Compare replay vs executed mid-course correction (and optional ceilings) on a campaign.

Parameters:: spec (EvaluateBatchControlPolicyInput)
Return type:: dict[str, Any]

process_improve.simulation.tools.get_simulation_tool_specs()[source]#

Return tool specs for all simulation tools registered in this module.

Return type:: list[dict]