Hard output constraints make a multi-fidelity reactor surrogate both data-efficient and safe to optimise against

Project identifier: SEPT-PCML Status: complete, all acceptance gates passed

Abstract

Surrogate models trained on expensive simulations are usually judged by prediction error, yet a surrogate that predicts well on held-out points can still be unusable for design because it violates conservation laws and because an optimiser will exploit whatever error it has. We asked whether enforcing stoichiometric conservation and non-negativity by construction, rather than by penalty, changes both the number of high-fidelity simulations required and the reliability of the resulting optimum. Using a non-isothermal Van de Vusse reactor as a frozen educational benchmark, we compared five model families across nested high-fidelity budgets of 32 to 512 samples, 30 seeds per configuration, and five evaluation regimes spanning interpolation, sparse regions, domain boundaries, near-thermal-limit states and declared out-of-distribution shells. The hard-constrained multi-fidelity model reduced the paired area under the learning curve by 45.2 percent relative to a black-box baseline (bootstrap 95 percent interval 43.9 to 47.6 percent). Constraint architecture separated the families more sharply than accuracy did. Hard-constrained models held the normalised stoichiometric residual at 4.4 x 10^-16, which is double-precision rounding, while the black-box and soft-penalty models reached residuals of 13.7 and 7.5 and emitted negative concentrations in up to 45.9 percent of predictions. The consequence appeared at the optimisation stage. Every surrogate optimum was re-evaluated in the high-fidelity simulator, and the black-box model predicted a yield of 1.0016 at its own optimum, which exceeds the physical maximum of unity, while the true verified yield was 0.6934. Verified regret ranged from 0.0745 for the soft-penalty model to 0.00087 for the hard-constrained multi-fidelity model. The advantage is not uniform: the single-fidelity hard model was more accurate on near-limit and out-of-distribution data, so we report a bounded rather than a general result. All values are numerical results for one frozen educational scenario and carry no claim of industrial reactor validity.

1. Introduction

Chemical reactor design depends on simulations that are accurate but slow. A surrogate model trained on a modest number of these simulations promises to make design search tractable, and this promise has driven wide adoption of machine-learning surrogates in process engineering. The standard evaluation reports a held-out error metric, and a low error is taken to indicate a usable model.

That evaluation misses two failures that matter for design. The first is physical admissibility. A surrogate that predicts species concentrations independently has no reason to conserve atoms, and nothing prevents it from returning a negative concentration. Such an output is not merely inaccurate, it is meaningless, and any downstream calculation built on it inherits the defect. The second failure is specific to optimisation. An optimiser searching a surrogate does not sample it uniformly. It is driven towards the regions where the surrogate is most optimistic, which are frequently the regions where the surrogate is most wrong. Average error over a test set says little about error at the point an optimiser will actually select.

Physics-informed machine learning addresses the first failure by adding residual penalties to the training loss. This approach reduces average violation, but a penalty is a preference rather than a guarantee, and the weight assigned to it trades against data fit. Whether soft penalties actually bring violations inside a tolerance, as opposed to merely shrinking their mean, is an empirical question that is often left unasked.

An alternative is to make the constraint part of the model's output parameterisation, so that no attainable set of network weights can produce a violating prediction. This removes the constraint from the loss entirely. The design cost is that the parameterisation must be derived for each conserved quantity, and it is not obvious in advance whether the restriction helps or hurts predictive accuracy, or whether any accuracy cost is repaid at the optimisation stage.

A separate consideration is the cost of the training data itself. When high-fidelity simulations are expensive, a multi-fidelity scheme that learns the discrepancy between a cheap approximate model and the expensive one can, in principle, reach a target accuracy with fewer expensive runs. Combining a multi-fidelity discrepancy model with a hard constraint layer raises a specific technical difficulty: adding an unconstrained correction to a constrained low-fidelity prediction destroys the constraint. The correction must therefore be learned in a latent space in which the constraint is preserved.

We evaluated these questions on a frozen benchmark built from a non-isothermal Van de Vusse reactor. The reaction network combines a series path and a parallel path, so that yield of the intermediate species is genuinely non-monotonic in residence time and temperature, which makes optimisation a real search rather than a boundary evaluation. Five model families were compared: an uncorrected low-fidelity mechanistic model, a black-box neural network, a soft-penalty network, a hard-constrained network, and a hard-constrained multi-fidelity model. Acceptance criteria, nested budgets, the split protocol and the optimisation objective were all frozen in a versioned specification before final evaluation began, so the comparison could not be tuned after the fact. These criteria were not registered on any external public registry.

Three questions organise the results. Does a hard constraint layer reduce the high-fidelity sample count needed to reach a given accuracy? Do soft penalties deliver the feasibility that hard parameterisation guarantees? And does either advantage survive the transition from prediction to verified optimisation?

2. Methods

2.1 Reactor model and scenario status

The high-fidelity model solves the non-isothermal Van de Vusse system in a constant-density time coordinate. The reaction network is

A -> B -> C          (series)
2A -> D              (parallel)

with rates r1 = k1(T) CA, r2 = k2(T) CB and r3 = k3(T) CA^2, and Arrhenius temperature dependence kj(T) = Aj exp(-Ej / RT). The species balances are

dCA/dtau = -r1 - 2 r3
dCB/dtau =  r1 - r2
dCC/dtau =  r2
dCD/dtau =  r3

and the energy balance, in one frozen sign convention, is

rho Cp dT/dtau = -sum_j (DeltaH_j r_j) - (UA/V)(T - Tc)

For this stoichiometry and a feed containing no products, the conserved inventory is

CA + CB + CC + 2 CD = CA0

This invariant is the constraint the hard models enforce exactly, and it holds only for the constant-density basis adopted here.

All kinetic and thermal parameter values are labelled EDUCATIONAL SCENARIO ASSUMPTION in configs/scenario_v1.yaml. They are internally consistent and support a well-posed numerical benchmark. They are not measured values for any real reactor, and the temperature bounds used in the output layer are numerical trust-domain limits rather than a safety envelope.

The low-fidelity model is the inherited isothermal Project H model, which is cheap and globally robust but structurally biased because it omits the thermal coupling entirely. Its regression suite of 136 tests is executed from an immutable snapshot as acceptance gate PML-G00.

2.2 Input and output definitions

Inputs are the vector x = [Tin, Tc, tau, CA0, UA_over_V]. The fundamental outputs are the four concentrations, the outlet temperature and the heat duty. Conversion, yield and selectivity are recomputed from the constrained fundamental outputs rather than predicted separately, which is what allows derived quantities to reconcile exactly. Selectivity is marked undefined below a frozen conversion tolerance instead of reporting a large and misleading ratio.

2.3 Dataset, splits and budgets

The dataset contains 1,792 matched low-fidelity and high-fidelity rows generated across five immutable strata: interior points, boundary faces and corners, near-thermal-limit points, sparse residence-time regions, and out-of-distribution shells reserved for stress testing. Generation recorded solver status, tolerances and balance residuals for every point, and produced zero failures.

Five disjoint test sets follow the strata. Preprocessing statistics are fit on training data only. A near-duplicate check across splits returned zero cross-split near duplicates, which is acceptance gate PML-G02. High-fidelity budgets are strictly nested at 32, 64, 128, 256 and 512 samples, so that each larger training set contains the smaller one and model comparisons at different budgets are paired rather than independent.

2.4 Model families

LF is the uncorrected low-fidelity mechanistic model, included as a bias reference.

BB is a black-box network trained on high-fidelity data with no physics term.

SOFT adds mass, energy, non-negativity and ODE residual penalties to the data loss. All penalty weights and scalings are recorded.

HARD applies the constraint in the output parameterisation. The network emits unconstrained logits z, and the composition is decoded as

p  = softmax(z)
CA = CA0 pA,  CB = CA0 pB,  CC = CA0 pC,  CD = (CA0 / 2) pD

The softmax guarantees non-negativity and sums to unity, so the conserved inventory CA + CB + CC + 2 CD = CA0 holds identically for any weights, up to floating-point rounding. Outlet temperature uses a bounded sigmoid transform into the declared trust domain, and heat duty is computed from the energy balance rather than predicted independently.

MF-HARD maps the low-fidelity prediction into the same logit space, adds the learned discrepancy there, and decodes through the identical hard layer:

latent_hi(x) = transform(low_fidelity(x), delta_theta(x))
y_hi_hat     = hard_projection(latent_hi, x)

Learning the correction in logit space rather than in concentration space is what preserves the constraint. An additive correction applied to concentrations directly would break the balance.

2.5 Training protocol and evaluation

Every family used identical splits, identical nested budgets, a fixed per-family tuning budget, and validation-only early stopping. The seed list was pre-declared. Thirty seeds were run per family and budget, giving 600 completed benchmark runs with zero failures. No architecture or hyperparameter was changed after final evaluation began.

Predictive accuracy is reported as a composite normalised RMSE across outputs. Physics metrics are the maximum and high-quantile normalised stoichiometric residual, the fraction of negative predictions, the energy-balance residual and the derived-output reconciliation residual. Uncertainty uses a 30-member ensemble with normalised split conformal calibration at nominal 0.90 coverage, calibrated on a set disjoint from both training and test data.

Data efficiency is measured as the paired reduction in area under the learning curve relative to BB, with a bootstrap 95 percent interval. The frozen criterion required a median improvement of at least 10 percent with a strictly positive lower bound.

2.6 Optimisation and verification

The frozen objective maximises yield of B subject to an outlet temperature ceiling of 400 K, an absolute heat duty limit of 8000 W, a minimum conversion of 0.7, and a requirement that selectivity be defined. Each surrogate was optimised inside the declared domain, and each resulting optimum was then re-evaluated in the high-fidelity simulator. Regret is the shortfall against a directly optimised high-fidelity reference of 0.7230. An optimum that proves infeasible on verification counts as a failure rather than as a low-error result.

3. Results

3.1 Hard constraints reduce the high-fidelity sample requirement

MF-HARD reached a given accuracy with fewer high-fidelity simulations than the black-box baseline at every nested budget (Figure 1a). The paired reduction in area under the learning curve was 45.2 percent, with a bootstrap 95 percent interval of 43.9 to 47.6 percent (Figure 1b). This clears the frozen threshold of 10 percent with a lower bound well away from zero, satisfying acceptance gate PML-G06.

The separation appears early. At 128 high-fidelity samples MF-HARD reached a median composite NRMSE of 0.0167, which the black-box model had not reached even at 512 samples, where it stood at 0.0199. On this benchmark the constrained multi-fidelity model reached with 128 simulations an accuracy the black-box model did not reach with four times as many.

Notably, the soft-penalty model tracked the black-box model closely throughout, with median IID errors of 0.0196 and 0.0199 respectively at the final budget. Adding physics residuals to the loss did not measurably improve sample efficiency here. The benefit came from the constraint architecture and the low-fidelity prior, not from the penalty.

Table 1. Median composite NRMSE by model family and evaluation regime at the final 512-sample budget (30 seeds).

Bold marks the best family in each column. MF-HARD is best on IID data, but HARD is best on near-limit and out-of-distribution data, and the two are close on sparse and boundary data.

The learning curves below show where that reduction comes from, and panel (b) shows the paired comparison the gate was frozen against.

3.2 Soft penalties reduce violations; hard parameterisation eliminates them

The constraint results separate the families far more sharply than the accuracy results do (Figure 2). Across all 600 runs, all five evaluation regimes and all budgets, the hard-constrained families held the maximum normalised stoichiometric residual at 4.34 x 10^-16 for HARD and 4.41 x 10^-16 for MF-HARD. These are double-precision rounding magnitudes, and they satisfy acceptance gate PML-G03.

Unconstrained families did not approach this. The black-box model reached a maximum normalised residual of 13.7 and the soft-penalty model reached 7.5. Both are of order unity or larger, meaning the predicted atom inventory departed from the feed inventory by more than the feed inventory itself in the worst case. The soft penalty roughly halved the worst-case violation relative to no penalty at all, which is a real effect, but it left the residual sixteen orders of magnitude above what the parameterisation achieves.

Negative concentrations show the same pattern in a more direct form. The black-box model returned at least one negative concentration in up to 44.9 percent of predictions within a run, and the soft-penalty model in up to 45.9 percent. The most negative single value observed was -6707 mol per cubic metre against a feed concentration of order 500. Both hard families returned exactly zero negative predictions, by construction.

Two further consistency results follow from the architecture. Because conversion, yield and selectivity are recomputed from constrained fundamentals, the maximum derived reconciliation residual was exactly 0 (gate PML-G04). Because heat duty is computed from the energy balance rather than predicted, the maximum hard-model energy residual was also exactly 0 (gate PML-G05).

This is the clearest result in the study. A soft penalty made violations smaller. It did not make them absent, and no penalty weight was found that did.

The two panels below separate the two behaviours: a penalty makes violations small, and a hard parameterisation makes them impossible.

3.3 The optimiser exploits whichever error the surrogate has

Unconstrained outputs show their consequences most starkly at the optimisation stage, where every surrogate optimum was verified in the high-fidelity simulator (Figure 3). All five optima were feasible on verification and none failed, satisfying acceptance gate PML-G08.

The black-box model predicted a yield of 1.0016 at its own optimum. A yield of B cannot exceed unity, because it is defined per mole of A fed, so the surrogate recommended an operating point on the basis of a physically impossible prediction. The verified yield at that point was 0.6934, an optimism gap of 0.308 yield units. The soft-penalty model showed the same pathology in weaker form, predicting 0.8779 against a verified 0.6485, a gap of 0.229.

Table 2. Optimisation results, all optima verified in the high-fidelity simulator.

High-fidelity reference optimum: 0.72300.

Ranking by regret does not match ranking by predictive accuracy. The soft-penalty model had almost the same IID error as the black-box model, 0.0196 against 0.0199, yet its verified regret was 2.5 times worse, 0.0745 against 0.0296. It was also worse than the uncorrected mechanistic model, whose regret was 0.0533. A model can be no less accurate on average and considerably more dangerous to optimise against, because what matters at the optimum is the error where the optimiser is pulled, not the error averaged over a test set.

MF-HARD achieved a verified regret of 0.000871, reaching a yield of 0.72213 against the reference 0.72300. Its optimism gap of 0.146 is not small. But the constrained output layer keeps the recommended point inside the physically admissible set, so the optimiser cannot escape into a region where the surrogate has invented performance no feasible state can deliver.

The figure below pairs each surrogate's own prediction with what the high-fidelity simulator returned at that same point, which is where the gap becomes visible.

3.4 The advantage is bounded, not general

No single family led everywhere (Figure 4a and Table 1). HARD was more accurate on near-limit data, 0.01698 against 0.03035, and on out-of-distribution data, 0.01465 against 0.02267. The low-fidelity prior helps where it is informative and hurts where its structural bias is largest, which is exactly at the thermal limits it omits by construction.

Uncertainty calibration was also regime-dependent (Figure 4b). Against a nominal 0.90 coverage, all families were adequately calibrated or conservative on IID and sparse data, where coverage ran from 0.917 to 0.992. Coverage degraded on boundary data, and MF-HARD was the worst affected at 0.844. Intervals were not recalibrated on test data, so this is reported as it stands rather than corrected. A user of this surrogate should treat its boundary intervals as optimistic.

These two results are why we report a bounded claim. Hard constraints eliminated violations at every point evaluated, across all 600 runs and all five regimes. The multi-fidelity accuracy advantage did not extend to every regime, and the interval calibration did not either.

Splitting the same models by evaluation regime shows the advantage changing hands, and the coverage panel shows where the intervals stop being trustworthy.

4. Discussion

The main result is that the choice between a soft penalty and a hard parameterisation is not a matter of degree. Across 600 runs and five evaluation regimes, hard parameterisation held conservation residuals at double-precision rounding while the soft penalty left worst-case residuals of order unity and negative concentrations in nearly half of its predictions. A penalty expresses a preference that competes with data fit. A parameterisation removes the possibility of violation from the model's range.

Why does this distinction matter beyond tidiness? The optimisation results answer that. A surrogate is often built to be searched, and search is adversarial towards surrogate error. The black-box model recommended an operating point on the basis of a predicted yield above the physical maximum, which is the signature of an optimiser finding and exploiting a region where the model has no physical grounding. The comparison between the soft-penalty model and the black-box model makes the point sharply, because their average accuracies were nearly identical while their verified regrets differed by a factor of 2.5. Held-out error did not predict optimisation reliability.

Scope matters for the data-efficiency result. A 45.2 percent paired reduction in area under the learning curve is substantial, and reaching at 128 simulations an accuracy the black-box model missed at 512 would matter when each simulation is expensive. But the mechanism is a useful low-fidelity prior combined with a constrained latent correction, and the benefit follows the prior. Where the isothermal low-fidelity model is structurally blind, at the thermal limits and outside the training distribution, the single-fidelity hard model was more accurate. A practitioner choosing between HARD and MF-HARD should decide on the regime they expect to operate in, not on the headline number.

Three limitations bound every claim above. First, all kinetic and thermal parameters are declared educational scenario assumptions. The benchmark is numerically well posed and internally consistent, and it is not a plant model; nothing here has been validated against experimental reactor data. Second, the temperature bounds in the output layer are numerical trust-domain limits, not safety limits, and should not be read as a safety envelope. Third, the comparison covers one reaction network, one constraint structure and one domain. Whether the ranking transfers to other conserved quantities or other reaction systems is untested.

The direction that would most strengthen the conclusion is a constraint structure where the conserved quantity is less convenient. The atom balance here maps onto a softmax with unusual neatness. Conservation laws that do not admit such a clean parameterisation are the real test of whether hard architectural constraints generalise as a design principle, and that test remains open.

5. Claims and non-claims

Supported by this work. All three claims below hold within the frozen educational Van de Vusse scenario and nowhere beyond it. The hard-constrained families satisfied stoichiometric conservation and non-negativity at double-precision tolerance everywhere evaluated. MF-HARD met the frozen paired data-efficiency criterion against the black-box baseline. All five surrogate optima were verified feasible in the high-fidelity simulator, with the regrets reported in Table 2.

Not supported by this work. Industrial reactor validity. Validation against plant data. Universal superiority of physics-constrained machine learning, which the near-limit and out-of-distribution results specifically contradict. Any physical safety envelope. Commercial reactor optimisation. The sufficiency of low prediction error as evidence of optimisation reliability, which Section 3.3 contradicts directly.