Scenario Quantile Predictive Control of Bioreactor Cultivation under Kinetic Parameter Uncertainty

2.1 Optimal Control of Cultivation Processes

Bioprocess trajectory optimization combines kinetic models, product-formation laws, and path constraints. Numerical methods for these problems are reviewed by Banga and co-workers [3]. The growth-associated and non-growth-associated terms of the Luedeking–Piret law provide a conventional product model [7]. Their use here specifies a simulator architecture; the nominal coefficients are study settings and are not estimates for a named organism.

2.2 Model Predictive, Robust, and Stochastic Control

Constrained MPC supplies a receding-horizon optimization framework [4], while stochastic MPC makes the probability model and constraint budget explicit [5]. The scenario approach yields probabilistic results for specified classes of sampled optimization problems, with assumptions and sample requirements tied to those results [8]. Those results do not automatically certify a nonconvex order-statistic program solved locally and reused in closed loop. The construction evaluated here therefore distinguishes its probability target from its empirical planning surrogate.

At the tested scale, tractability depends on horizon length, number of scenarios, and nonlinear-solver behavior. Warm starts reduce repeated optimization effort, but a comparison of one nominal model with ten scenarios is not a controlled scaling experiment.

2.3 Parameter Uncertainty and Estimation in Bioprocesses

Kinetic estimation and adaptive bioreactor control have an established literature [6]. Bayesian inference provides a distribution over uncertain parameters [9], and sequential Monte Carlo supplies one possible update mechanism [10]. A learned distribution requires an explicit likelihood, assay design, identifiability assessment, and predictive checks. None of these inference steps is exercised in the reported campaign, which samples a known simulation law.

2.4 Synthesis and Research Gap

The present study connects uncertain kinetics to predictive substrate control through an empirical quantile constraint. Its contribution is the paired comparison and quantile ablation under a specified simulator, rather than a claim that posterior-aware bioprocess control has no precedent. Table 1 distinguishes the assumptions of established control families from the empirical status of this implementation.

3.1 Research Questions and Hypotheses

The study asks how scenario-based kinetic uncertainty affects the productivity–violation trade-off and whether the resulting optimization is inexpensive at the tested horizon. The questions below distinguish a design target from an observed result.

RQ1: Can an empirical scenario quantile provide a tunable substrate-risk mechanism in the reported simulator, and how does its observed violation frequency relate to the nominal probability target?

RQ2: How do titer and substrate excursions differ among the six controllers under paired parameter draws and process-noise seeds?

RQ3: What solve-time budget is observed for the stated horizon and scenario count?

The titer hypothesis compares the quantile controller with the fixed-margin controller, while the risk hypothesis examines whether a higher scenario order statistic reduces observed violations. The computation hypothesis concerns the measured budget at the executed scale.

3.2 System Model and Notation

The general fed-batch model has biomass, substrate, and product concentrations together with culture volume, x = (X, S, P, V) and a clean feed of rate F and substrate concentration S_in. Eq. (1) includes the dilution of biomass and product caused by increasing volume. The illustrative inhibited growth law in Eq. (2) uses θ = (μ_max, K_S, K_I, K_P, Y, a, b). The executed comparison instead uses the constant-volume system of Eq. (6), a six-parameter law without substrate inhibition. Table 2 separates the physical quantities from the controller settings.

3.3 Formal Problem Statement

At each control cycle, the ideal stochastic program of Eq. (4) balances an expected objective with marginal state constraints of Eq. (3). Its probability must cover the parameter and disturbance laws used to describe the plant. A joint path requirement would instead constrain the event that all inequalities hold simultaneously and would need a separate batch or horizon risk budget. The reported controller approximates the marginal target using deterministic parameter scenarios; it does not solve the ideal program exactly.

dX F = (μ − )X, dt V dS μX F = − + (Sin − S), dt Y V dP F = (aμ + b)X − P, dt V dV = F.

dt S KP

μ(S, P) = μmax S2 · KP + P. (2)

KS + S + K

Prθ,w{gj(xk+ℓ) ≤ 0} ≥ 1 − εj, for each l and j. (3) minF Eθ,w[J(x, F; θ, w)], subject to (1), (3), 0 ≤ F ≤ Fmax, V ≤ Vmax.

Π(θ | y) ∝ p(y | θ)p0(θ). (5)

4.1 Posterior Construction for Kinetic Parameters

A posterior-informed implementation would combine a kinetic likelihood with prior information and update it when assays become available [9], [10]. Such a layer would need parameter identifiability and posterior-predictive checks before its uncertainty could support a control decision. The executed study uses a specified lognormal law directly and contains no assay-driven inference. This distribution is therefore a known simulation input, not an estimated posterior.

4.2 Scenario-Based Chance-Constrained MPC

At every control cycle, the implemented controller propagates ten deterministic parameter scenarios over the prediction horizon. The scenario set is drawn once per run and reused, as described in Section 5.1. For each horizon step, the ninth ordered substrate prediction must not exceed the ceiling; a quadratic mean-exceedance penalty also enters the objective. This is an empirical nearest-rank 90% planning rule. It permits one of ten scenarios to exceed the ceiling and provides no established 90% population coverage guarantee. The convex scenario results of [8] require assumptions and sample bounds that have not been established for this locally solved nonlinear quantile program.

The scenario order statistic is applied separately at each predicted step. A pointwise scenario constraint does not certify the simultaneous event over the full horizon or batch, and process noise is absent from the planning scenarios. The empirical step and run violation summaries in Section 5 measure these distinct effects. They are diagnostics of the chosen design, not substitutes for a closed-loop probabilistic proof.

4.3 Comparison Controllers

The nominal MPC benchmark uses one nominal kinetic vector. Hard-scenario MPC enforces the substrate ceiling for every sampled parameter vector. The fixed-margin benchmark, labelled tube-lite in the reported tables and figure, is nominal MPC with a 15% substrate back-off. It contains no propagated tube, affine disturbance feedback, or min-max optimization. Its zero observed violations must therefore be interpreted as a result of this campaign, not as recursive feasibility or worst-case safety.

4.4 Numerical Implementation

The executed controllers share explicit-Euler prediction and warm-started SLSQP single shooting through SciPy [15]. These are locally solved nonlinear programs [11]; global optimality and recursive feasibility are not established. Previous control blocks initialize the next solve. Reproducible implementation requires the exact objective and terminal surrogate, control-block parameterization, stopping tolerances, unsuccessful-solve policy, and timing protocol. The reported aggregate solve times appear in Section 5.2.

4.5 Constraint Handling, Complexity, and Failure Modes

The general design distinguishes actuator and vessel limits from recoverable state limits. Its executed instantiation constrains only the dilution-rate actuator and substrate ceiling. Increasing the number of scenario rollouts adds prediction work, but nonlinear-solver iterations also affect total cost. Prior misspecification, poor identifiability, and structural model error are validation concerns for a future learned-parameter implementation; this campaign does not quantify them.

The actuator constraint is 0 ≤ F ≤ F_max in the general fed-batch model and 0 ≤ u ≤ 0.5 h⁻¹ in the executed constant-volume model. The volume ceiling V ≤ V_max belongs to the general model and is inactive in the

5.1 Experimental Setup

The reported campaign instantiates the control design on a constant-volume, three-state bioreactor surrogate with x = (X, S, P) and manipulated dilution rate 0 ≤ u ≤ 0.5 h⁻¹. The dilution losses in Eq. (6) correspond to a well-mixed continuous-culture model with equal inflow and outflow. They do not describe a changing-volume fed-batch vessel. Growth follows Monod kinetics [1] with multiplicative product inhibition, and product formation uses the Luedeking–Piret structure [7]:

dX = (μ − u)X, dt dS μX = u(Sin − S) − , dt Y

dP (6)

= (βμ + γ)X − uP, dt

μ = μmax .

KS + S P KP The stated discretization is explicit Euler with Δt = 0.2 h for 60 steps, giving T = 12.0 h. The feed concentration is S_in = 20.0 g/L and x₀ = (0.2, 2.0, 0.0) g/L. Nominal coefficients are μ_max = 0.4 h⁻¹, K_S = 0.8 g/L, Y = 0.5 g/g, β = 1.5 g/g, γ = 0.02 g/(g h), and K_P = 6.0 g/L. The substrate ceiling is 4.0 g/L with a nominal 0.9 marginal probability target. Volume, product, and terminal biomass guards are inactive. Because Eq. (6) has no Andrews term [2], the ceiling does not mark a modeled transition into substrate inhibition.

The six parameters follow the study-selected law θ = θ_nom exp(0.2z), z ~ N(0,I₆), independently by component. The nominal vector is the median, not the mean; each component has coefficient of variation √(exp(0.04) − 1), approximately 20.2%. Evaluation crosses 30 parameter draws with five process-noise seeds. Within each of the 150 cells, controllers use the same true parameter vector and noise realization. Additive per-step Gaussian noise has standard deviations (0.01, 0.01, 0.02) g/L, followed by truncation at zero. The two scenario controllers share ten planning parameter draws sampled once per run. A violation is S > 4.0 g/L at a sampled step. The campaign reports a total runtime of 168.2 s on a desktop Central Processing Unit (CPU).

OL control solves the nominal single-shooting program once (reported solve time 114 ms) with four 15-step constant-control segments, producing u = (0.036, 0.019, 0.250, 0.150) h⁻¹. PI control [16] tracks 3.9 g/L with conditional integration; K_p = 0.2 and K_i = 0.005 were selected from a stated 5 × 5 nominal-plant search. Nominal MPC uses H = 10 steps and two control blocks. Hard-scenario MPC (SMPC) enforces the maximum substrate prediction over ten scenarios. The proposed quantile MPC (CCMPC) uses the ninth order statistic and a quadratic mean squared-exceedance penalty with ρ = 2.0. Fixed-margin MPC (TUBE in the tables) imposes predicted S ≤ 3.4 g/L. All MPC variants use the stated horizon-end-product and remaining-production surrogate; its exact formula is required to reproduce the comparison. The Chance-Constrained Model Predictive Control (CCMPC)– Stochastic Model Predictive Control (SMPC) comparison changes both the constraint and the penalty, so it is not a pure quantile ablation.

The reported outcomes are final titer P(T), productivity P(T)/T, sampled-step violation fraction, peak substrate, and solve time. Tables give mean ± standard deviation, rather than confidence intervals, across the 150 paired cells; lower-decile titer is the mean of the lowest 10% of runs. The reported software versions are Python 3.12.15, NumPy 2.1.3 [17], SciPy 1.14.1 [15], and Matplotlib 3.9.2 [18]. Figures are reported as generated from result files at 300 dpi.

The implemented study isolates the control layer under a known parameter law. It uses full-state feedback and a shared Euler discretization for plant and planner, with parametric mismatch and additive noise. It does not test posterior estimation, structural mismatch, min-max or certified tube control, changing-volume fed-batch operation, or inhibition- and biomass-specific guards. The quantile ablation uses 15 draws crossed with three seeds, a 45-cell subset.

5.2 Main Comparison (RQ1, RQ2)

Table 3 reports the six-method comparison over the 150 paired runs; Fig. 1 shows substrate trajectories against the chance-constraint bound for the three highest-titer controllers on a representative draw, and Fig. 2 places every individual run in the violation-titer plane. Three regimes are immediately visible.

The nominal OL profile has a reported predicted peak of 4.000 g/L and nominal titer of 4.507 g/L. Applied to the uncertain systems, it gives 2.788 ± 1.130 g/L, or 61.9% of that nominal prediction, and violates the ceiling in 49.1% of sampled steps and all 150 runs. The reported mean peak substrate is 12.22 g/L and the maximum is 14.841 g/L, about 3.7 times the ceiling. The early low feed and mid-course rise expose the sensitivity of an unchanged profile to kinetic mismatch. These results show overfeeding relative to the operational guard; they do not quantify biological damage from substrate inhibition, which is absent from Eq. (6).

PI and fixed-margin MPC have no observed violations in this campaign, with reported peaks of 3.908 and 3.483 g/L. Their mean titers are 4.798 ± 1.516 and 4.913 ± 1.506 g/L. Relative to CCMPC, their differences are −0.251 and −0.136 g/L, corresponding to discounts of approximately 5.0% and 2.7% using CCMPC as denominator. PI has a lower-decile titer of 2.334 g/L but a highest draw-mean titer of 8.467 g/L. Thus the classical loop is competitive within the simulated conditions. Neither zero observed violations nor a fixed planning margin constitutes a safety guarantee.

The mean titers of Nonlinear Model Predictive Control (NMPC), SMPC, and CCMPC differ by at most 0.015 g/L, while pooled step violation fractions are 0.065, 0.032, and 0.049 and runs with any violation are 33.3%, 16.7%, and 25.3%, respectively. NMPC is a point-parameter controller; the other two propagate parameter scenarios. SMPC has the higher reported mean and lower violation frequency than CCMPC, so these data do not establish a performance advantage for the quantile form. Mean solve times are 3.10, 3.90, 4.40, and 3.98 ms for NMPC, SMPC, CCMPC, and fixed-margin MPC. The 4.40-ms value is about 164,000 times shorter than the 720-s sampling interval. This is a measured implementation budget, not a real-time worst-case bound; solver failure and timing variability require separate analysis.

Lines: Seed-mean paths of the three highest-titer controllers; shaded bands: Min-max envelope across the five process-noise seeds; dashed line: S_max (Fig. 1).

Fig. 1 shows similar trajectories for the three MPC controllers on draw 21. The paths approach the ceiling, and the seed envelope includes shallow excursions above it. Their similarity is consistent with the small titer differences in Table 3; the aggregate violation summaries address behavior beyond this illustrative draw.

Each method has 150 simulated cells. The tube-lite label denotes fixed-margin nominal MPC (Fig. 2).

Fig. 2 places PI and fixed-margin MPC on the observed zero-violation axis, the three MPC variants at higher mean titer with nonzero excursions, and OL control at a substantially larger violation frequency. This comparison is an empirical trade-off across the tested designs, rather than a certified productivity–safety frontier.

5.3 Ablation: The Quantile Dial

Table 4 compares nearest-rank levels 80%, 90%, and 95% on the stated 45-cell subset. For ten scenarios, the selected ranks are ceil (10q), namely eight, nine, and ten. The 95% rule therefore uses the scenario maximum; its difference from SMPC includes the additional penalty. On the same subset, SMPC has reported titer 5.470 ± 1.610 g/L, violation fraction 0.031, eight violating runs, and lower-decile titer 3.235 g/L.

Observed violation fractions decrease from 0.079 to 0.051 to 0.029 as the selected rank increases, and violating-run fractions decrease from 33.3% to 26.7% to 17.8%. Mean titers are 5.434, 5.482, and 5.421 g/L, so the response is not monotone. The reported results include run-level Wilcoxon values W = 0, p = 6.5 × 10⁻⁴ for 80% versus 95% violations; W = 435, p = 0.358 and W = 423, p = 0.292 for 90% versus 80% and 95% titer; and W = 3, p = 2.8 × 10⁻¹³ for 80% versus 95% titer, with median difference 0.006 g/L and 44 of 45 pairs aligned. The last comparison is statistically different under the reported calculation, although its mean gap is only 0.013 g/L. The three seeds within each of 15 draws limit independent inference. These results support an empirical quantile response, not a universally optimal 90% setting.

The constraint controls an order statistic of scenario substrate predictions at each step. The quadratic exceedance penalty acts on squared positive violations and is not Conditional Value-at-Risk (CVaR), whose optimization construction is different [19]. The 95% row and SMPC have similar violation fractions, 0.029 and 0.031, and a reported titer difference of −0.049 g/L. Without a separate penalty-coefficient ablation and solver diagnostics, the results do not identify the causal contribution or gradient strength of the penalty.

5.4 Robustness across Parameter Draws

Table 5 reports parameter-draw means over the five noise seeds. CCMPC draw-mean titer spans 1.775–8.269 g/L. The reported analysis identifies draw 29 with μ_max = 0.243 h⁻¹ and draw 13 with β = 2.725 g/g, and reports titer correlations of +0.84 with μ_max, −0.53 with K_S, +0.48 with β, −0.03 with Y, and −0.12 with K_P. These are descriptive associations across 30 draws and do not isolate parameter effects or prove that the same uncertainty controls titer and constraint risk. The substrate guard does not protect final titer against adverse kinetics.

Among runs with any violation, reported mean violating-step fractions are 0.194, 0.189, and 0.193 for NMPC, SMPC, and CCMPC, corresponding to approximately 11–12 sampled steps or 2.3 h. CCMPC violations arise within 15 draws; SMPC has 19 fully violation-free draws. The reported largest MPC peak is 4.058 g/L, approximately 1.5% above the ceiling, compared with 14.841 g/L for OL control. These durations and magnitudes help distinguish excursions, but their mechanisms require trajectory-level examination rather than inference from aggregates alone.

The pooled step satisfaction fractions are 93.5% for NMPC, 95.1% for CCMPC, and 96.8% for SMPC. A pooled average over times and runs does not establish Pr{S(t) ≤ S_max} ≥ 0.9 at every time. Fully violation-free runs account for 66.7%, 74.7%, and 83.3%, respectively, and therefore do not meet a 90% batch-level empirical target. On the ablation subset, the 95% order statistic raises violation-free runs from 73.3% to 82.2%. Fixed-margin MPC has zero observed violations, but none of these observations certifies a joint guarantee for future batches.

5.5 Statistical Analysis

The analysis uses Wilcoxon signed-rank comparisons [20], a Friedman screen [21], and Holm correction across five baseline comparisons [22]. Table 6 reports those calculations on the 150 paired run cells. The parameter-draw unit is shared by five noise seeds, however, so the cells are nested rather than 150 independent biological draws. The reported p-values are exploratory; confirmatory inference requires a draw-level analysis or a justified clustered procedure.

Table 6 reports a Friedman statistic of 359.1, and CCMPC paired mean differences of +2.261, +0.251, +0.005, −0.010, and +0.136 g/L against OL, PI, NMPC, SMPC, and fixed-margin MPC. Its primary CCMPC–SMPC calculation gives W = 5348 and p = 0.555; this absence of rejection does not establish equivalence. The analysis additionally reports a post hoc violation comparison, W = 172.5 and p = 8.2 × 10⁻⁴, with 38 versus 25 violating runs. The independent-draw limitation prevents treating them as definitive significance or superiority evidence.

5.6 Interpretation of the Reported Comparison

The comparison supports three descriptive findings. An unchanged nominal OL profile performs poorly under the stated uncertainty. Feedback controllers reduce substrate excursions substantially, and the scenario representations have similar mean titers but different observed violation frequencies. Within the quantile ablation, a higher selected rank reduces observed violations with small mean-titer changes. The measured optimization cost is modest at this scale. The 26.1% difference between SMPC and NMPC timing does not establish a sublinear scenario-count scaling law, because their constraints and solver paths also differ.

SMPC is ahead of CCMPC in the reported mean titer, violation fraction, and lower-decile titer. Fixed-margin MPC and PI have no observed violations and achieve approximately 97.3% and 95.0% of CCMPC mean titer; PI wins 46 of 150 paired cells. These observations make baseline selection consequential. A practitioner would need application-specific valuations of excursions, computation, and product output before choosing a controller; the simulation does not designate a universally preferable policy.

The known-law scenarios remove estimation error from the experiment, but they do not provide a mathematical optimistic bound on safety under a learned law. Common plant and planner discretization also excludes structural error and integration mismatch. The evidence concerns one simulator, one dispersion scale, one horizon, and the stated paired draws. Its positive result is an empirical mechanism for adjusting observed substrate risk; certified guarantees, posterior learning, and organism-specific validation remain separate tasks.

6.1 Limitations

The study uses a continuous-culture surrogate rather than a volume-explicit fed-batch model. It omits substrate inhibition, dissolved-oxygen and pH interactions, sensor faults, contamination, and estimator delay. Thus the substrate ceiling has an operational role, and final concentration is not a complete measure of harvested mass or economic productivity. Euler step-size convergence and unsuccessful-solve behavior also need verification. The six-component lognormal law is a simulation assumption, not a measured biological variability model.

6.2 Alternative Viewpoints and Responses

The empirical comparison illustrates when simple policies can be competitive. PI provides direct substrate regulation, and a fixed planning margin reduces observed excursions without scenario propagation. A posterior-aware extension could use estimation tools [6], [9], [10], but credible risk statements would require the resulting distribution to be validated and the control and disturbance assumptions to be stated. Scenario-count, prior-misspecification, structural-mismatch, and penalty-coefficient studies would separate the mechanisms that the present comparison changes together.