Paper deep dive
A tale of perfect fit and phantom optima: how data-driven models can fail in real-time optimization
Prithvi Dake, Rahul Bindlish, James B. Rawlings
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/26/2026, 4:56:48 AM
Summary
This paper investigates the reliability of data-driven models for Real-Time Optimization (RTO) in chemical processes. Using a vinyl acetate monomer benchmark process, the authors demonstrate that both structured hybrid models (combining mass balances with neural network closures) and fully data-driven neural ODE models can achieve accurate fits to historical plant data while producing significantly different and suboptimal economic optima compared to the actual plant. The study highlights that these 'phantom optima' arise not only from model structure but also from the stochastic nature of the training optimizer, which can drift weights away from the true optimum even with noise-free data. The authors propose a decision-oriented benchmark requiring models to recover the plant optimum, rather than just fitting data, before plant application.
Entities (7)
Relation Signals (5)
Good Predictive Fit → doesnotguarantee → Reliable Economic Performance
confidence 95% · These results demonstrate that a good predictive fit of all available measurements does not guarantee reliable economic performance.
Vinyl Acetate Monomer Process → servesas → Decision-Oriented Benchmark
confidence 92% · A data-driven model for RTO should at least be required to recover the optimum on a decision-oriented benchmark like the one developed here
Structured Hybrid Model → produces → Phantom Optima
confidence 90% · Yet their economic optima differ substantially from that of the plant. Where the plant returns a single optimum on multistart search, the trained models return many phantom optima.
Neural Ordinary Differential Equation → produces → Phantom Optima
confidence 90% · Yet their economic optima differ substantially from that of the plant. Where the plant returns a single optimum on multistart search, the trained models return many phantom optima.
Stochastic Gradient Training → causes → Suboptimal RTO Solutions
confidence 88% · Even with noise-free data and initialization at weights that recover the plant optimum, stochastic gradient training can drift to weights that yield substantially worse RTO solutions.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Real-time optimization (RTO) relies on process models to locate economically optimal operating conditions. Because developing first-principles models requires significant process knowledge, data-driven alternatives are increasingly attractive. Modern machine-learning models can fit historical plant data accurately and often pass standard validation tests. Whether such models can be trusted for economic optimization, however, remains unclear. We investigate this question using a vinyl acetate monomer benchmark process with a unique, well-conditioned economic optimum. We train a structured hybrid model that combines known mass balances and thermodynamics with a neural-network closure for unknown kinetics, and a fully data-driven neural ordinary differential equation (ODE) model. Both models reproduce plant measurements accurately and exhibit little variation in predictions across random initializations. Yet their economic optima differ substantially from that of the plant. Where the plant returns a single optimum on multistart search, the trained models return many phantom optima. We further show that the training optimizer alone can be yet another source of error. Even with noise-free data and initialization at weights that recover the plant optimum, stochastic gradient training can drift to weights that yield substantially worse RTO solutions. The identified model is thus an artifact of the training optimizer as well as the data. These results demonstrate that a good predictive fit of all available measurements does not guarantee reliable economic performance. A data-driven model for RTO should at least be required to recover the optimum on a decision-oriented benchmark like the one developed here before being considered for plant testing and application.
Tags
Links
- Source: https://arxiv.org/abs/2608.23885v1
- Canonical: https://arxiv.org/abs/2608.23885v1
Trouble viewing inline? Open PDF directly →
Full Text
84,383 characters extracted from source content.
Expand or collapse full text
A tale of perfect fit and phantom optima: how data-driven models can fail in real-time optimization Journal: Computers and Chemical Engineering Prithvi Dake Email: prithvidake@ucsb.edu Corresponding author: Corresponding Author Address: Department of Chemical Engineering, University of California, Santa Barbara, CA 93106, United States Rahul Bindlish Email: rbindlish@dow.com Address: Dow Chemical Company, TX, United States James B. Rawlings Email: jbraw@ucsb.edu Address: Department of Chemical Engineering, University of California, Santa Barbara, CA 93106, United States Abstract Real-time optimization (RTO) relies on process models to locate economically optimal operating conditions. Because developing first-principles models requires significant process knowledge, data-driven alternatives are increasingly attractive. Modern machine-learning models can fit historical plant data accurately and often pass standard validation tests. Whether such models can be trusted for economic optimization, however, remains unclear. We investigate this question using a vinyl acetate monomer benchmark process with a unique, well-conditioned economic optimum. We train a structured hybrid model that combines known mass balances and thermodynamics with a neural-network closure for unknown kinetics, and a fully data-driven neural ordinary differential equation (ODE) model. Both models reproduce plant measurements accurately and exhibit little variation in predictions across random initializations. Yet their economic optima differ substantially from that of the plant. Where the plant returns a single optimum on multistart search, the trained models return many phantom optima. We further show that the training optimizer alone can be yet another source of error. Even with noise-free data and initialization at weights that recover the plant optimum, stochastic gradient training can drift to weights that yield substantially worse RTO solutions. The identified model is thus an artifact of the training optimizer as well as the data. These results demonstrate that a good predictive fit of all available measurements does not guarantee reliable economic performance. A data-driven model for RTO should at least be required to recover the optimum on a decision-oriented benchmark like the one developed here before being considered for plant testing and application. Keywords: Real-time optimization , Physics-constrained machine learning , Neural ordinary differential equations , Parameter estimation , Decision-oriented benchmark 1 Introduction Growing global competition and sustainability demands have made RTO increasingly important to the chemical process industries (Naysmith and Douglas, 1995; Darby et al., 2011; Câmara et al., 2016). RTO operates on a timescale of hours to days, or is triggered when a price change or disturbance is detected. It supplies the profitable steady-state setpoints to lower-level regulatory layers (Cutler and Perry, 1983). Traditional RTO strategies often use a rigorous nonlinear steady-state process model that is updated using plant measurements after each optimization iteration (Seborg et al., 2017, pg. 351). Specifically, the model parameters are reconciled against the most recent steady-state plant measurements, a step termed model adaptation or data reconciliation (Crowe, 1996). The updated model is then used to solve a steady-state economic optimization problem. This two-step approach (Roberts and Williams, 1981) is a common way of practicing RTO. While the plant tracks the RTO setpoint, it is in transient operation, and one must detect the new steady state before reconciling with measurements. Thus, the approach is limited by steady-state wait-times and fragile detection algorithms that require careful tuning (Darby et al., 2011). Recent developments in the two-step RTO approach adapt the steady-state model using transient plant measurements (Krishnamoorthy et al., 2018), which eliminates both the steady-state wait-times and detection. The RTO problem can also be made dynamic, e.g., economic MPC (Rawlings et al., 2012). However, when the disturbances vary slowly relative to the plant’s settling time, steady-state RTO is preferred because it is cheaper to solve. Besides model-based methods, model-free methods also exist. RTO has been carried out using direct-search methods that usually mimic gradient-free schemes like Nelder-Mead or evolutionary strategies (Box, 1957). Without gradient information, the number of plant evaluations grows steeply with the decision space for these methods. A second family of model-free methods can be understood as a feedback problem in the form of extremum seeking or self-optimizing type control (Skogestad, 2000). If we treat the gradient of profit with respect to the degrees of freedom at steady-state as a controlled variable and hold it at zero, the plant self-optimizes to a stationary point of the RTO problem. Most of the feedback methods therefore focus on gradient estimation, e.g., using dither signals. The feedback approach does not distinguish among a maximum, minimum, or saddle point, and leads to prohibitively slow convergence and a difficult dither tuning problem in practice (Krishnamoorthy and Skogestad, 2022). For these reasons, model-free methods are not considered further. A third method, modifier adaptation, combines the two: it converges in few plant iterations, like a model-based method, yet provably reaches the plant’s optimum, like the feedback-based schemes (Bonvin and Pannocchia, 2024). Crucially, it recovers the plant optimum despite model mismatch, provided the model already has the correct curvature at the plant optimum (a positive semidefinite Hessian in the case of a minimum), so that the plant’s gradient corrections suffice (Chachuat et al., 2009; Papasavvas et al., 2019). Yet this curvature is not guaranteed even for a model with no mismatch, as we show in this paper. The accurate estimation of plant gradients remains another challenge. We therefore focus on model-based RTO using transient measurements. RSDChemical Plant(A)(B)Plant Data Figure 1: (A) Historical plant datasets are used to train (B) a structured model and a black-box ML/AI closure model, each deployed in real-time optimization (RTO). (C) Near-identical fits map to very different optima, the worst-case structured model loses close to 30% profit and the black-box model over 50%. A reliable data-driven model must recover the optimum, not be judged solely on the validation dataset. 1.1 Modeling For RTO, practitioners usually start with first-principles models (FPMs), built on conservation laws and thermodynamics. The model’s structure leaves few free parameters, so it extrapolates well beyond the operating data and remains interpretable. These are exactly the properties an economic optimizer relies upon (Bindlish, 2025). But the cost to build and maintain such models is large since FPMs demand more process knowledge than is often available. An FPM is only as good as its assumed structure. When that structure is incomplete, e.g., when the reaction kinetics or fluid dynamics are not known from first principles, the model violates the adequacy condition of Forbes et al. (1994): no choice of parameters satisfies the plant’s optimality conditions, and its economic optimum is biased. Alternatives include purely data-driven or black-box models that learn plant dynamics from input–output data with little structure imposed. Such models use neural networks (NNs) that are universal function approximators (Hornik et al., 1989). Neural ordinary differential equations (Chen et al., 2018) are a representative modern instance of such black-box models. They need little process knowledge, but they do not preserve conservation laws and degrade rapidly outside the training set. Violating conservation laws leads to infeasible predictions, which is usually the biggest concern in chemical engineering applications. Augmenting the first-principles balances with a data-driven closure for the unknown parts in FPMs has therefore been a long-standing pursuit in chemical engineering (see the reviews by Sansana et al., 2021; Schweidtmann et al., 2024). The resulting modeling techniques have been recently surveyed under the umbrella of physics-constrained machine learning (PCML) by Mukherjee and Zavala (2026). PCML strategies are commonly grouped by how the physics is enforced. In soft-constrained models, the known physical laws enter the loss as a penalty, as in physics-informed neural networks (Raissi et al., 2019; Sholokhov et al., 2023). In contrast, hard-constrained models enforce steady-state mass and energy balances exactly (Golder et al., 2025; Constante et al., 2026). Another example includes the simultaneous projection method, which is seen as a promising step toward using NNs to approximate parts of DAE models (Lueg et al., 2025). Finally, in hybrid (or structured) models, a data-driven term is wired into a partially mechanistic transient model and solved with shooting methods (Kumar and Rawlings, 2023; Thompson et al., 2026). The physics-based dynamic balances are retained, while a data-learned term stands in for the part that is not known from first principles, typically the reaction kinetics or fluid dynamics. We use both hybrid and black-box models in this work. Recent perspectives highlight the absence of decision-oriented benchmarks as central open problems in PCML (Pantelides et al., 2025; Mukherjee and Zavala, 2026). At the same time, a growing literature argues that reliable models for process systems can be built largely from historical plant data using machine learning (ML) and artificial intelligence (AI) (Sansana et al., 2021; Schweidtmann et al., 2024). We take up both threads and ask: for a data-driven model built with the best modeling paradigms and software available to us, carrying as much physical structure as the process permits, and clearing the usual validation checks, does the optimum it yields under RTO match the plant’s? A graphical abstract of the work is shown in fig. 1. Our contributions are as follows: (i) Evidence that near-identical fits, tight even across random initializations, still yield sharply disparate economic optima, so standard validation and uncertainty quantification do not certify reliable RTO performance. (i) A controlled demonstration that the stochastic training optimizer alone can cause divergence from the optimal weights. Given correct model structure, noise-free data, and initialization at weights at which the model recovers the plant optimum, stochastic gradient training still drifts to a model whose RTO solution is significantly suboptimal. (i) A decision-oriented benchmark, the vinyl acetate process with a well-conditioned economic optimum, on which a data-driven RTO model must recover the plant optimum rather than merely fit plant data. Outline The paper is organized as follows: Section 2 introduces the vinyl acetate process used as the case study, sets up its economic optimization, and describes the structured (hybrid) and black-box (neural ODE) models. Section 3 gives the training data, scaling, initialization, and software details. Section 4 reports the results, comparing the parametric, structured, and black-box models on how well they recover the plant economic optimum. Section 5 concludes. Further details on the flowsheet are provided in Appendices A and B. The supplementary material reports the hyperparameter tuning and the training optimizer comparison for every case studied in this paper. 2 Industrial case-study: Vinyl acetate process The vinyl acetate monomer (VAc) production process has been used as a benchmark problem for plantwide control study due to its challenging and coupled process network (Luyben and Tyréus, 1998; Chen et al., 2003). We model the process using a detailed reactor model coupled with a simplified separator system, which includes the absorber, flash drum and distillation column. Although the separator system is simplified, the nonlinear reactor model provides rich insight into the process for steady-state economic optimization (Ward et al., 2004). 2.1 Process description The VAc process involves seven chemical species: ethylene (C2H4C X_ [t]2H X_ [t]4) with inert ethane (C2H6C X_ [t]2H X_ [t]6), oxygen (O2O X_ [t]2) and acetic acid (CH3COOHCH X_ [t]3COOH) as fresh feeds; vinyl acetate (C4H6O2C X_ [t]4H X_ [t]6O X_ [t]2) is the product; water (H2OH X_ [t]2O) and carbon dioxide (CO2CO X_ [t]2) are byproducts. The following reactions take place inside the reactor: Wedevelopasimplifiedflowsheetasshownin.Intheprocess,thechemicalspecies(ingaseousphase)comingfromthemixer(stream2)arefedtothereactor.Thegaseouseffluentfromthereactor(stream3)ispartlycondensedtoseparateintogasandliquidintheflashdrum.Thegases(stream5)thenpassthroughaCO 2 absorber,andstream6isthenpartlyrecycledbacktothemixerbythesplitterI(stream7).Theliquidstream4isthensenttoadistillationcolumn.Theoverheadproducts(mainlyvinylacetateandwater)arewithdrawn(stream10)whilethebottomproducts(mainlyaceticacidandwater)arerecycledbacktothemixerbysplitterII(stream13).Thus,theflowsheetconstitutesaconventionalreactor-separatorsystemwithtworecycles.Notethatthereactionsinsection 2.1arenotelementarybutratheranoverallrepresentationoftheactualreactionnetwork.Wealsoassumealumpedmodelofthecatalyst,suchthatallthemasstransferresistancesarerepresentedbythefollowingratelaws: r 1 = k 1 × e - E 1 / R T × c Ee 1.1 c O 1.1 c A r 2 = k 2 × e - E 2 / R T × c Ee c O wherec_j, Taretheconcentrationofspeciesaretheconcentrationofspeciesjandthetemperatureinsidethereactor.Thenetproductionratesaregivenas, (1) Equation 1 eq. Equation eqs. Equations 1 = [ R Ee R Ea R A R V R W R O R C ] [ - 1 - 1 0 0 - 1 0 1 0 1 2 - ( / 1 2 ) - 3 0 2 ] [ r 1 r 2 ] Wemodelthereactorasasinglecontinuousstirred−tankreactor(CSTR)operatingisothermallywithnegligiblepressuredrop.Thereactorismaintainedatthesametemperatureasinletstream2in,whosetemperatureisasetpointwespecify.Weassumethatthethermodynamicequationofstateisvalidevenwhenthereactorisnotatequilibrium(rawlings:ekerdt:2020, pg. 133).Assuminganidealgasmixture,thecomponentbalancesforspeciesandthetemperatureinsidethereactor.Thenetproductionratesaregivenas, equation bmatrixR_Ee\\ R_Ea\\ R_A\\ R_V\\ R_W\\ R_O\\ R_C bmatrix= bmatrix-1&-1\\ 0&0\\ -1&0\\ 1&0\\ 1&2\\ -(1/2)&-3\\ 0&2 bmatrix bmatrixr_1\\ r_2 bmatrix equationWemodelthereactorasasinglecontinuousstirred-tankreactor(CSTR)operatingisothermallywithnegligiblepressuredrop.Thereactorismaintainedatthesametemperatureasinletstream2in @crefcreftype~refnumfig:block_measurement,whosetemperatureisasetpointwespecify.Weassumethatthethermodynamicequationofstateisvalidevenwhenthereactorisnotatequilibrium [citep](\@@bibrefAuthorsPhrase1Yearrawlings:ekerdt:2020\@@citephrase, , pg. 133).Assuminganidealgasmixture,thecomponentbalancesforspeciesjinthereactorare:ϵVRdcjdt=Qfcjf−Qcj+(1−ϵ)VRρcRjQ=Qf+RTP(1−ϵ)VRρc∑jRj Figure 2 Figure 2 fig. Figure figs. Figures 2 Figure 2 7 $ $ $ $ $ $ Gas recycle α Reactor N 1 Ee N 1 A N 1 O Flash Distillation column Absorber Mixer Splitter I 1 2 8 3 4 10 11 9 6 5 12 Splitter I 13 Liquid recycle β c Ee c Ea c A c W c V c O c C TC 2Simplified process flowsheet for the vinyl acetate (VAc) monomer plant consisting of 7 units and 13 streams. (Ee: Ethylene, Ea: Ethane, A: Acetic acid, O: Oxygen, W: Water, V: Vinyl acetate, C: Carbon dioxide). Outlet concentrations from the reactor are measured subject to step inputs in the decision variables (fresh feed rates, recycle fractions and reactor temperature), that are shown by the control valves and the controller icon. Figure 2Simplified process flowsheet for the vinyl acetate (VAc) monomer plant consisting of 7 units and 13 streams. (Ee: Ethylene, Ea: Ethane, A: Acetic acid, O: Oxygen, W: Water, V: Vinyl acetate, C: Carbon dioxide). Outlet concentrations from the reactor are measured subject to step inputs in the decision variables (fresh feed rates, recycle fractions and reactor temperature), that are shown by the control valves and the controller icon. Wemakethefollowingassumptionstosimplifytheflowsheet:ethane,ethylene,oxygenandcarbondioxidearenon−condensableandthus,wegetperfectcondensationofonlyaceticacid,waterandvinylacetateintheflashdrum.Wealsoassumethattheseparationinthedistillationcolumnissharp.Thus,stream10carriesallthevinylacetateandwaterproducedinthereactorwhilestream8containsonlyaceticacid.Hence,neithervinylacetatenorwatergetsrecycledtothemixer.TheCO2absorberisassumedtoremoveonlycarbondioxidefromstream5.Finally,forunitoperationsotherthanthereactor,onlymaterialbalancesaresolved;energybalancesarenotconsidered.Letinthereactorare: aligned ε V_R dc_jdt&=Q_fc_jf-Qc_j+(1-ε)V_R _cR_j\\ Q&=Q_f+ RTP(1-ε)V_R _c _jR_j aligned figure* 15826875 picture (0.24575,0.61422) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l7 tabular (0.0425,0.45273) [rgb]1,0,0 @makebox@ (0.0,0.0)[lt] tabular[t]l\$ tabular (0.23142,0.6526) [rgb]1,0,0 @makebox@ (0.0,0.0)[lt] tabular[t]l\$ tabular (0.25539,0.06829) [rgb]1,0,0 @makebox@ (0.0,0.0)[lt] tabular[t]l\$ tabular (0.90215,0.21603) [rgb]0,1,0 @makebox@ (0.0,0.0)[lt] tabular[t]l\$ tabular (0.0,0.0) [width=627.37706pt,page=1]Economics-fig.pdf (0.90125,0.64845) [rgb]1,0,0 @makebox@ (0.0,0.0)[lt] tabular[t]l\$ tabular (0.0,0.0) [width=627.37706pt,page=2]Economics-fig.pdf (0.8977,0.07611) [rgb]1,0,0 @makebox@ (0.0,0.0)[lt] tabular[t]l\$ tabular (0.0,0.0) [width=627.37706pt,page=3]Economics-fig.pdf (0.49014,0.63702) [rgb]1,0,0 @makebox@ (0.0,0.0)[lt] tabular[t]lGas recycle $α$ tabular (0.0,0.0) [width=627.37706pt,page=4]Economics-fig.pdf (0.43993,0.32389) @makebox@ (0.0,0.0)[t] tabular[t]cReactor tabular (0.01151,0.36245) @makebox@ (0.0,0.0)[t] tabular[t]c$N^1_Ee$ tabular (0.00411,0.32824) @makebox@ (0.0,0.0)[t] tabular[t]c$N^1_A$ tabular (0.00575,0.29334) @makebox@ (0.0,0.0)[t] tabular[t]c$N^1_O$ tabular (0.0,0.0) [width=627.37706pt,page=5]Economics-fig.pdf (0.75209,0.32839) @makebox@ (0.0,0.0)[t] tabular[t]cFlash tabular (0.0,0.0) [width=627.37706pt,page=6]Economics-fig.pdf (0.75574,0.17297) @makebox@ (0.0,0.0)[t] tabular[t]cDistillation\\ column tabular (0.0,0.0) [width=627.37706pt,page=7]Economics-fig.pdf (0.75311,0.47619) @makebox@ (0.0,0.0)[t] tabular[t]cAbsorber tabular (0.0,0.0) [width=627.37706pt,page=8]Economics-fig.pdf (0.18598,0.32474) @makebox@ (0.0,0.0)[t] tabular[t]cMixer tabular (0.0,0.0) [width=627.37706pt,page=9]Economics-fig.pdf (0.75218,0.5968) @makebox@ (0.0,0.0)[t] tabular[t]cSplitter I tabular (0.0,0.0) [width=627.37706pt,page=10]Economics-fig.pdf (0.04737,0.40788) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l1 tabular (0.27595,0.34114) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l2 tabular (0.75875,0.09444) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l8 tabular (0.61692,0.3415) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l3 tabular (0.75775,0.24366) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l4 tabular (0.89997,0.17517) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l10 tabular (0.90771,0.49513) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l11 tabular (0.91296,0.61515) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l9 tabular (0.75595,0.53042) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l6 tabular (0.75563,0.40727) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l5 tabular (0.0,0.0) [width=627.37706pt,page=11]Economics-fig.pdf (0.89885,0.04317) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l12 tabular (0.75574,0.02418) @makebox@ (0.0,0.0)[t] tabular[t]cSplitter I tabular (0.0,0.0) [width=627.37706pt,page=12]Economics-fig.pdf (0.23649,0.03837) [rgb]0,0,1 @makebox@ (0.0,0.0)[lt] tabular[t]l13 tabular (0.47588,0.05626) [rgb]1,0,0 @makebox@ (0.0,0.0)[lt] tabular[t]lLiquid recycle $β$ tabular (0.0,0.0) [width=627.37706pt,page=13]Economics-fig.pdf (0.62421,0.31137) @makebox@ (0.0,0.0)[t] tabular[t]c$c_Ee$ tabular (0.62711,0.28415) @makebox@ (0.0,0.0)[t] tabular[t]c$c_Ea$ tabular (0.6204,0.25158) @makebox@ (0.0,0.0)[t] tabular[t]c$c_A$ tabular (0.62371,0.21938) @makebox@ (0.0,0.0)[t] tabular[t]c$c_W$ tabular (0.61941,0.18927) @makebox@ (0.0,0.0)[t] tabular[t]c$c_V$ tabular (0.62083,0.15583) @makebox@ (0.0,0.0)[t] tabular[t]c$c_O$ tabular (0.61756,0.12353) @makebox@ (0.0,0.0)[t] tabular[t]c$c_C$ tabular (0.0,0.0) [width=627.37706pt,page=14]Economics-fig.pdf (0.26347,0.20448) [rgb]0.102,0.102,0.102 @makebox@ (0.0,0.0)[lt] tabular[t]lTC tabular picture \@@toccaption @tag[ ]2Simplified process flowsheet for the vinyl acetate (VAc) monomer plant consisting of 7 units and 13 streams. (Ee: Ethylene, Ea: Ethane, A: Acetic acid, O: Oxygen, W: Water, V: Vinyl acetate, C: Carbon dioxide). Outlet concentrations from the reactor are measured subject to step inputs in the decision variables (fresh feed rates, recycle fractions and reactor temperature), that are shown by the control valves and the controller icon.\@@caption @tag[: ]Figure 2Simplified process flowsheet for the vinyl acetate (VAc) monomer plant consisting of 7 units and 13 streams. (Ee: Ethylene, Ea: Ethane, A: Acetic acid, O: Oxygen, W: Water, V: Vinyl acetate, C: Carbon dioxide). Outlet concentrations from the reactor are measured subject to step inputs in the decision variables (fresh feed rates, recycle fractions and reactor temperature), that are shown by the control valves and the controller icon. \@add@centering figure*Wemakethefollowingassumptionstosimplifytheflowsheet:ethane,ethylene,oxygenandcarbondioxidearenon-condensableandthus,wegetperfectcondensationofonlyaceticacid,waterandvinylacetateintheflashdrum.Wealsoassumethattheseparationinthedistillationcolumnissharp.Thus,stream10carriesallthevinylacetateandwaterproducedinthereactorwhilestream8containsonlyaceticacid.Hence,neithervinylacetatenorwatergetsrecycledtothemixer.TheCO X_ [t]2absorberisassumedtoremoveonlycarbondioxidefromstream5.Finally,forunitoperationsotherthanthereactor,onlymaterialbalancesaresolved;energybalancesarenotconsidered. Let N_j^k bethemolarflowofchemicalspeciesbethemolarflowofchemicalspeciesj ∈J ≔ Ee, Ea, A, V, W, O, C instreaminstreamk ∈K ≔1, 2, …, 13 .Hence,themolarflowcomingoutofthereactorcanbewrittenas.Hence,themolarflowcomingoutofthereactorcanbewrittenasN_j^3 = Q c_j.Sinceonlythereactorcontributestodynamicsintheflowsheet,allotherprocessunitsareassumedtooperateatquasisteady−state.Thus,thematerialbalanceforthemixeratanytimecanbeconciselyrepresentedas:N2=N1+N7+N13=N1+diag(ζ)⋅N3ζ=(α,α,β,0,0,α,α(1−μ))TQf=∑jNj2RT/Pwhere.Sinceonlythereactorcontributestodynamicsintheflowsheet,allotherprocessunitsareassumedtooperateatquasisteady-state.Thus,thematerialbalanceforthemixeratanytimecanbeconciselyrepresentedas: aligned N^2&=N^1+N^7+N^13=N^1+diag(ζ)· N^3\\ ζ&=(α,α,β,0,0,α,α(1-μ))^T\\ Q_f&= _jN_j^2RT/P alignedwhereα≔N_j^7 / N_j^6, β≔N_j^13/N_j^8aretherecyclefractionscorrespondingtostreams7and13respectively,whilearetherecyclefractionscorrespondingtostreams7and13respectively,whileμ≔N_C^11/N_C^5isthefractionofmolarflowofcarbondioxideremovedbytheCO2absorber.TheconcentrationsatthereactorinletcanbecalculatedusingisthefractionofmolarflowofcarbondioxideremovedbytheCO X_ [t]2absorber.Theconcentrationsatthereactorinletcanbecalculatedusingc_jf = N_j^2 / Q_f.Sections 2.1,,2.1and2.1describethetransientbehavioroftheflowsheetformingadifferentialalgebraicequation(DAE)system.WecansimplifytheDAEintoanODEbysolvingsection 2.1for. @crefcreftypepluralcap~refnumeq:rxnrate, @crefrefnumeq:netprodrate, @crefrefnumeq:mass_baland @nobreakspace @crefrefnumeq:mixer_baldescribethetransientbehavioroftheflowsheetformingadifferentialalgebraicequation(DAE)system.WecansimplifytheDAEintoanODEbysolving @crefcreftype~refnumeq:mixer_balforQ_fintermsofintermsofc_jonly.Forthesteady−stateoftheflowsheet,weneedonlytosettheRHStermsoftheODEequationsequaltozero.Notethatthefeedsonly.Forthesteady-stateoftheflowsheet,weneedonlytosettheRHStermsoftheODEequationsequaltozero.NotethatthefeedsN_W^1, N_V^1, N_C^1aresettozerowhiletheinertfeedentersasaresettozerowhiletheinertfeedentersasN_Ea^1 = λN_Ee^1wherewhereλisasmallconstantnumber.ConstructingtheremainingstreamsintheflowsheetisdescribedinAppendix A.isasmallconstantnumber.Constructingtheremainingstreamsintheflowsheetisdescribedin @crefcreftypecap~refnumapp:process_units. 2.2 Economic optimization A candidate objective function that takes into account the steady-state economics of the flowsheet is defined as follows, ℓ(xs,us,p)≔ (x_s,u_s;p) p1NV10−p2∑jNj7−p3∑jNj13 p_1N_V^10-p_2 _jN_j^7-p_3 _jN_j^13 −∑jp4jNj1−p5∑jNj9−p6∑jNj12 - _jp_4jN_j^1-p_5 _jN_j^9-p_6 _jN_j^12 (42) which is used to formulate the following RTO problem: maxxs,usℓ(xs,us,p) _x_s,u_s (x_s,u_s;p) s.t.f(xs,us)= s.t. 9.24994ptf(x_s,u_s)= 0 0 g(xs,us)= g(x_s,u_s)= NV,nom10−NV10≤0 N^10_V,nom-N_V^10≤ 0 u¯≤us≤u¯ u≤ u_s≤ u (43) where x are the states which are the species concentrations c in the reactor and c≔(cEe,cEa,cA,cV,cW,cO,cC)Tc (c_Ee,c_Ea,c_A,c_V,c_W,c_O,c_C)^T. The decision variables (or optimization degrees of freedom) u are the feed rates, recycle fractions and the reactor inlet temperature, u≔[NEe1,NO1,NA1,α,β,T]Tu [N^1_Ee,N^1_O,N^1_A,α,β,T]^T The price vector p≔(p1,p2,p3,p4jT,p5,p6)Tp (p_1,p_2,p_3,p_4j^T,p_5,p_6)^T characterizes the profit function ℓ(⋅) (·). The RHS of the ODE formed by sections 2.1, , 2.1 and 2.1 is represented by f(⋅)f(·). The constraint g(⋅)g(·) ensures that the nominal production rate of vinyl acetate is met. We set NV,nom10=0N^10_V,nom=0 to make the problem unconstrained. We also fix the liquid recycle β=0.6β=0.6 to remove the flat direction it forms with NA1N^1_A. The vector p and the bounds on the degrees of freedom u¯,u¯ u, u are tabulated in Appendix B. The subscript s denotes the steady-state values of the states and the degrees of freedom. The meaning of the terms in the objective function is as follows: • p1p_1 represents the revenue from the vinyl acetate production NV10N_V^10. While the remaining costs represent the overall operating expenditure (opex). • p2,p3p_2,p_3 represent energy processing costs for gas and liquid recycle streams N7,N13N^7,N^13, while p5,p6p_5,p_6 represent the cost to process the purge streams N9,N12N^9,N^12. • p4jp_4j represent the material costs associated with each species j in the fresh feed N1N^1. If we supply excessively high feed rates, the residence time in the reactor and consequently the yield of vinyl acetate will be low. On the other hand, very low feed rates will lead to low production rates. A high gas recycle can lead to a better yield but incurs a large energy cost. Thus, an unconstrained optimum should exist in terms of the degrees of freedom for the above objective function. Solving eq. 43, we obtain us⋆=[15.54mol/s,11.14mol/s,17.15mol/s,0.57,0.6,427.1K]Tu_s =[15.54\ mols,11.14\ mols,17.15\ mols,0.57,0.6,427.1\ K]^T with ℓ⋆=292.87$/sNV10=13.13mol/s =292.87\ \$/s 18.49988ptN_V^10=13.13\ mols We analyze the reduced Hessian of the profit function with respect to the degrees of freedom. Because u contains variables with different units, we rescale them using u~ u with u~=u/su u=u/s_u, where the scale factors sus_u are tabulated in Appendix B, so that the condition number reflects curvature rather than the choice of units. The eigenvalues of the scaled reduced Hessian ∂2ℓ∂u~s2 ∂^2 ∂ u_s^2 are all negative, confirming a local maximum, and its condition number is 415, so the optimum is locally well-conditioned with no flat directions. We performed a multi-start search from different initial guesses for usu_s and xsx_s. It revealed no other local optima within u¯≤us≤u¯ u≤ u_s≤ u. Together with the negative definite scaled reduced Hessian, this establishes a strict, well-conditioned local maximum with no competing optimum detected in the feasible region. At the computed optimum, profit is 7.4%7.4\% of revenue, a margin typical of commodity chemicals, so the optimum reflects a genuine economic trade-off rather than a tiny difference between two large numbers. We retain the decision variable scaling for the models developed in the upcoming sections. Problem eq. 43 is solved using the IPOPT solver provided by CasADi (Wächter and Biegler, 2006; Andersson et al., 2019). 2.3 Structured models We consider sections 2.1, , 2.1 and 2.1 as the ground truth for the flowsheet. In reality, however, some parts of the plant model are not known precisely. The reactor is where the chemical transformation takes place, and the kinetics are usually not known from first principles. This is almost always the case in industrial practice (Bindlish, 2025). We assume that the separation factors for the simple separators are known. Thus, we choose to parameterize only the rate laws using feedforward neural networks as follows: r1 r_1 =rNN,1(cEe,cO,cA,T,θrNN,1) =r_N,1(c_Ee,c_O,c_A,T; _r_N,1) r2 r_2 =rNN,2(cEe,cO,T,θrNN,2) =r_N,2(c_Ee,c_O,T; _r_N,2) (44) where θrNN,1,θrNN,2 _r_N,1, _r_N,2 represent the weights and biases of the neural networks. These rate laws replace section 2.1 in the plant model and are trained as follows. If we represent the ODE formed by sections 2.1 and 2.1 by x˙=f(x,u) x=f(x,u), then we can write the following optimization problem to train the neural-network parameters ΘN≔(θrNN,1,θrNN,2) _N ( _r_N,1, _r_N,2) as, ΘN⋆ _N =argminΘN(1NtrNt∑i=1Ntr∑k=0Nt|yi(k)−y^i(k)σy|2)1/2 = _ _N ( 1N_trN_t _i=1^N_tr _k=0^N_t | y_i(k)- y_i(k) _y |^2 )^1/2 s.t.x^˙i s.t. 9.24994pt x_i =f(x^i,ui,ΘN) =f( x_i,u_i; _N) y^i y_i =h(x^i,ui) =h( x_i,u_i) x^i(0) x_i(0) =f0(yi(0)) =f_0(y_i(0)) (45) where x,ux,u have been defined in section 2.2 and y=xy=x, thus h(x,u)=Ih(x,u)=I. NtN_t represents the number of time-steps in a trajectory while NtrN_tr represents the number of trajectories in the training dataset. Training over a prediction horizon (NtN_t) reduces the compounding of one-step errors during open-loop rollout (Kumar and Rawlings, 2023). The symbol σy _y represents the standard deviation of the measurements y over the training dataset, while the hat represents model predictions. For full-state measurement, we simplify the training by initializing each trajectory at its first measurement rather than estimating the initial conditions, i.e., f0(yi(0))=yi(0)f_0(y_i(0))=y_i(0). Since the model shown in eq. 45 enforces both mass conservation and the equation of state, we call it a structured model. As written, problem eq. 45 is the most realistic formulation in an industrial setting, where only concentration measurements are available. However, we consider one more way to train the networks in eq. 44. Suppose that we have access to the rate measurements [r1,r2]T[r_1,r_2]^T over the same training dataset used for eq. 45, thus we can train the neural networks simply as, ΘN⋆ _N =argminΘN(1NtrNt∑i=1Ntr∑k=0Nt|ri(k)−r^i(k)σr|2)1/2 = _ _N ( 1N_trN_t _i=1^N_tr _k=0^N_t | r_i(k)- r_i(k) _r |^2 )^1/2 (46) where σr _r is the standard deviation of the rate measurements over the training dataset. We do not know the neural-network parameters that reproduce the plant kinetics exactly. A rate measurement reflects the kinetic rate itself, rather than its integrated effect on concentration through the ODE in eq. 45. Hence eq. 46 gives the closest estimate of the closure we can obtain from data. Training on rates therefore acts as both a reference and a diagnostic. It tests whether the network has the capacity to represent the kinetics and whether the resulting model recovers the plant’s RTO optimum. If it does, we treat the rate-trained network as the best attainable surrogate for the true closure and judge the concentration-trained models against it. If it does not, the shortfall lies with the dataset rather than the training, since no measurement is more informative about the kinetics than the rates themselves, and the dataset must be enlarged. Thus using eqs. 45 and 46 we consider three cases, defined together with all other cases studied in table 1: (a) rmeas: Problem eq. 46 with access to rate measurements [r1,r2]T[r_1,r_2]^T. (b) cmeas_rinit: Problem eq. 45 initialized with the N parameters learned from rmeas. (c) cmeas: Problem eq. 45 with access to only concentration measurements c. 2.4 Black-box models We also train a complete black-box model for the reactor as, θfNN⋆ θ _f_N =argminθfNN(1NtrNt∑i=1Ntr∑k=0Nt|yi(k)−y^i(k)σy|2)1/2 = _ _f_N ( 1N_trN_t _i=1^N_tr _k=0^N_t | y_i(k)- y_i(k) _y |^2 )^1/2 s.t.x^˙i s.t. 9.24994pt x_i =fNN(x^i,ui,θfNN) =f_N( x_i,u_i; _f_N) y^i y_i =h(x^i,ui) =h( x_i,u_i) x^i(0) x_i(0) =f0(yi(0)) =f_0(y_i(0)) (47) where the RHS f of the ODE is approximated using a neural network. As in the previous section, we use y=xy=x and set f0(yi(0))=yi(0)f_0(y_i(0))=y_i(0). By construction, eq. 47 preserves neither individual species mass conservation nor the equation of state inside the reactor; it learns them from data instead. Hence, we call it a black-box model. To use it for RTO, we need to calculate the flowsheet streams. For that, we need to know the volumetric flow Q coming out of the reactor. From section 2.1, by applying an overall mass-balance across the flowsheet at steady state we get, Q=∑jMjNj1∑jMj(1−ζj)cj Q= _jM_jN^1_j _jM_j(1- _j)c_j (48) where MjM_j is the molecular weight of the species j∈Jj∈ J. Thus, using eqs. 47 and 48 and Nj3=QcjN^3_j=Qc_j we can construct the flowsheet and use the black-box model for economic optimization. 3 Training and implementation details 3.1 Training dataset As shown in , we record the species concentrations at the reactor outlet while perturbing the decision variables with pseudo-random binary signals under a zero-order hold. The decision variables are the feed rates NEe1,NO1,NA1N^1_Ee,N^1_O,N^1_A, the recycle fractions α,βα,β and the reactor inlet temperature T. We sample 200 random inputs and hold out 30% as a validation set. The inputs u, except temperature, are varied by ±30%± 30\% around their optimal values us⋆u_s , while temperature is varied by ±10%± 10\%. These limits coincide exactly with the box constraints of the RTO problem eq. 43, i.e., u¯=0.7us⋆ u=0.7\,u_s and u¯=1.3us⋆ u=1.3\,u_s for all the feed rates and recycle fraction α, and u¯=0.9T⋆ u=0.9\,T , u¯=1.1T⋆ u=1.1\,T for the temperature, tabulated in Appendix B. Although β is perturbed ±30%± 30\% to excite the data, it is held fixed at 0.60.6 in the RTO problem eq. 43 (i.e., u¯=u¯ u= u for β). The training dataset therefore spans exactly the feasible box over which we optimize. These perturbations are more aggressive than industrial model identification typically permits, so the training data here is richer than a practitioner could likely obtain. A subset of the noise-free training dataset is shown in fig. 3, and it brackets the measurements at the steady-state optimum xs,usx_s,u_s. Thus, we test the models purely on interpolation and not extrapolation. We ensure that the training data has enough steady-state information by holding the decision variables constant for a fixed period before injecting the next random step, so that the model learns the steady-state behavior of the flowsheet (Kumar and Rawlings, 2023). Figure 3: A snapshot of training data from the plant model with no measurement noise. The measurements are shown in red and the decision variables in black. The dotted grey lines represent xs⋆x_s and us⋆u_s from section 2.2. 3.2 Software details Both the structured and black-box models are built and trained using the JAX open-source Python library developed by Google DeepMind (Bradbury et al., 2018). We solve the ODEs using Diffrax (Kidger, 2021) with an adaptive step-size solver of the ESDIRK5(4) type. The neural networks are constructed and initialized using Flax (Heek et al., 2024). The hyperparameters in this paper are the network width and depth, the learning rate, the batch size, and the epoch budget. We fixed the width and depth because the training error was already low at this architecture, so additional capacity was unnecessary. We also fixed the batch size and tuned only the learning rate and the epoch budget. We first swept the learning rate and kept the value with the lowest validation loss, then increased the epoch budget until the validation loss plateaued. The tuning procedure for each model we construct is detailed in to of the supplementary material. We run ADAM (Kingma and Ba, 2014) for the first half of the epoch budget and then switch to L-BFGS (Liu and Nocedal, 1989) for the second half for faster convergence. A comparison of using only ADAM versus ADAM followed by L-BFGS is also shown in the supplementary material. Both the optimizers are obtained from Optax provided by JAX. 3.3 Scaling and initialization We use hyperbolic tangent as the activation function for the neural networks. Since both structured and black-box models are used inside a nonlinear program (NLP) eq. 43, tanh ensures differentiability of the model outputs with respect to the inputs and parameters. For tanh, both the inputs and outputs of the neural networks should be scaled so that the network operates on (1)O(1) quantities. The inputs are normalized by their mean and standard deviation over the entire training dataset. For the structured model, we deliberately separate the shape of the rate function from its scale: the network learns the standardized shape, while a trainable location μ^r μ_r and scale σ^r σ_r restore the physical magnitude. For example, r1r_1 is scaled as r1=μ^r1+σ^r1×rNN,1(⋅), r_1= μ_r_1+ σ_r_1× r_N,1(·), so that rNN,1(⋅)r_N,1(·) need only output an (1)O(1) standardized rate, while μ^r1 μ_r_1 and σ^r1 σ_r_1 carry its mean and spread. This scaling keeps the raw network output well-scaled and the back-propagated gradients well-conditioned; without it, the network would have to emit values spanning the raw rate magnitudes directly, which readily causes vanishing or exploding gradients, especially for stiff systems (Fronk and Petzold, 2025). To further stabilize training, we apply gradient clipping to the network parameters ΘN _N, but leave the scale factors μ^ri,σ^ri μ_r_i, σ_r_i unclipped. The weights are thereby protected from large, destabilizing updates, while μ^ri μ_r_i and σ^ri σ_r_i remain free to take the larger steps needed to reach the correct rate magnitude. We apply clipping only during the stochastic ADAM phase. The deterministic L-BFGS phase needs none, since its strong Wolfe line search already bounds each step size. The factors μ^ri,σ^ri μ_r_i, σ_r_i must therefore be initialized with care. For problem eq. 45, we initialize at the mean and standard deviation of [r1,r2]T[r_1,r_2]^T computed from the plant model over the training dataset, so that the network begins by learning only the standardized shape. In practice, the plant rates are unavailable. One can instead use eq. 47 to get the RHS and back-calculate approximate rates using section 2.1 to set these initial estimates as shown by Peng et al. (2025), though that is not the focus of this work. For problem eq. 46, by contrast, we do not need these scale factors, since they are already available from the assumed rate measurements. For the black-box model, we use a similar scaling strategy. The inputs to the neural network fNNf_N are normalized by their mean and standard deviation over the training dataset. The outputs are scaled by trainable location and scale factors μ^fNN,σ^fNN μ_f_N, σ_f_N that restore the physical magnitude while the network learns the standardized shape of the RHS of the ODE. As done for the structured model, we initialize μ^fNN μ_f_N and σ^fNN σ_f_N at the mean and standard deviation of the true RHS x˙ x computed from the plant model over the training dataset. Table 1: Summary of the six data-driven cases. nln_l is the sensor-noise level in eq. 50. A random N initialization uses the LeCun normal initialization popular for tanh-based neural networks (LeCun et al., 2002, pg. 20). Model Case Training problem and measurements N initialization Noise nln_l Structured rmeas eq. 46, rates random 00 cmeas_rinit eq. 45, concentrations rmeas solution 00 cmeas eq. 45, concentrations random 00 cmeas_noise eq. 45, concentrations random 10−210^-2 Black-box bbox eq. 47, concentrations random 00 bbox_noise eq. 47, concentrations random 10−210^-2 4 Results and discussion We demonstrated in section 2.2 that the flowsheet defined by sections 2.1, , 2.1 and 2.1 has a well-conditioned optimum. As a control experiment, we first verify that the plant rate structure is identifiable from the data the plant generated, and that the identified model recovers the plant optimum. This rules out uninformative data or an ill-posed problem as a trivial explanation for any failure of the structured or black-box models. Whether the optimum remains recoverable without the correct structure is precisely the question these experiments address. 4.1 Control experiment using parametric model From section 2.1, we note that the plant rate laws are power law models. We rescale and parameterize them as follows: r1= r_1= k~1×e−E~1/R(1/T−1/Tm) k_1× e^- E_1/R(1/T-1/T_m) ×(cEe/cEe,m)a1(cO/cO,m)b1(cA/cA,m)c1 ×(c_Ee/c_Ee,m)^a_1(c_O/c_O,m)^b_1(c_A/c_A,m)^c_1 r2= r_2= k~2×e−E~2/R(1/T−1/Tm) k_2× e^- E_2/R(1/T-1/T_m) ×(cEe/cEe,m)a2(cO/cO,m)b2 ×(c_Ee/c_Ee,m)^a_2(c_O/c_O,m)^b_2 (49) where Tm,cj,mT_m,c_j,m are the mean temperature and concentrations over the subset used to fit the parametric model eq. 49. This subset is roughly 0.25%0.25\% of the entire training-and-validation dataset used for the neural networks, shown in fig. 3. From it we estimate the parameters [logk~1,logE~1,a1,b1,c1,logk~2,logE~2,a2,b2][ k_1, E_1,a_1,b_1,c_1, k_2, E_2,a_2,b_2] along with their confidence intervals. To avoid biasing the estimator, we initialize far from the true plant parameters. We also corrupt each measurement with multiplicative Gaussian sensor noise, c~=c(1+ηnl)=c+ηnlc, c=c\,(1+η\,n_l)=c+η\,n_l\,c, (50) where η∼(0,1)η (0,1) is drawn independently for each measurement and nl∈0, 10−3, 10−2n_l∈\0,\,10^-3,\,10^-2\ is the relative error level (0, 0.1%, 1%). The parameter estimation problem is solved using Paresto (Dake et al., 2024) which provides joint confidence intervals for the parameters. Paresto solves a similar problem as eq. 45 using IPOPT, with r1,r2r_1,r_2 replaced by eq. 49 and ΘN _N by parameters in eq. 49. The results are shown in table 2. We see that the plant rate structure is identifiable from the training data, as indicated by the tight confidence intervals for the parameters even in the presence of noise. Each of the fits reports a single optimum even on multi-start RTO. We measure this recovery with the following metric, henceforth called the profit loss, Δℓ=ℓu^s−ℓ⋆ℓ⋆×100% = _ u_s- × 100\% (51) where ℓ⋆ is the optimal plant profit and ℓu^s _ u_s is the profit obtained by evaluating the plant model at the estimated operating point u^s u_s from an identified model. Note that Δℓ≤0% ≤ 0\% since ℓ⋆ is the global optimum. We see from table 2 that the profit loss is negligible at all noise levels, which indicates that the identified parametric model is able to recover the plant optimum. Table 2: Parameter estimation results and profit loss for the control experiment with parametric models for different levels of measurement noise. Parameter nl=0n_l=0 nl=10−3n_l=10^-3 nl=10−2n_l=10^-2 logk~1 k_1 −4.836-4.836 ± 2.7×10−112.7× 10^-11 −4.840-4.840 ± 4.4×10−34.4× 10^-3 −4.876-4.876 ± 3.9×10−23.9× 10^-2 logE~1 E_1 10.6110.61 ± 4.3×10−114.3× 10^-11 10.6010.60 ± 7.1×10−37.1× 10^-3 10.5510.55 ± 6.7×10−26.7× 10^-2 a1a_1 1.1001.100 ± 5.2×10−115.2× 10^-11 1.0921.092 ± 8.6×10−38.6× 10^-3 1.0211.021 ± 7.7×10−27.7× 10^-2 b1b_1 1.1001.100 ± 5.5×10−115.5× 10^-11 1.0911.091 ± 9.0×10−39.0× 10^-3 1.0141.014 ± 8.1×10−28.1× 10^-2 c1c_1 1.0001.000 ± 4.4×10−114.4× 10^-11 0.9930.993 ± 7.2×10−37.2× 10^-3 0.9340.934 ± 6.5×10−26.5× 10^-2 logk~2 k_2 −7.992-7.992 ± 1.1×10−111.1× 10^-11 −7.993-7.993 ± 1.8×10−31.8× 10^-3 −7.994-7.994 ± 1.8×10−21.8× 10^-2 logE~2 E_2 10.8710.87 ± 1.8×10−111.8× 10^-11 10.8710.87 ± 3.0×10−33.0× 10^-3 10.8710.87 ± 3.0×10−23.0× 10^-2 a2a_2 1.0001.000 ± 2.4×10−112.4× 10^-11 0.99960.9996 ± 4.0×10−34.0× 10^-3 0.9960.996 ± 4.0×10−24.0× 10^-2 b2b_2 1.0001.000 ± 4.0×10−114.0× 10^-11 1.0001.000 ± 6.7×10−36.7× 10^-3 0.99970.9997 ± 6.7×10−26.7× 10^-2 Δℓ <−0.001%<-0.001\% −0.019%-0.019\% −1.28%-1.28\% These results can be used as a benchmark to compare the performance of the structured and black-box models. 4.2 Training and validation of the data-driven models (a) Structured models Because neural networks are over-parameterized and their training loss is non-convex, each random initialization of the weights and biases can settle in a different basin and recover a distinct rate function. For each of the three cases rmeas, cmeas_rinit and cmeas, we therefore train an ensemble of 20 networks, each from an independent random initialization. Note that we randomly initialize only the N weights and biases and not the scaling factors. We train one additional case called cmeas_noise, which is similar to cmeas but uses noisy measurements eq. 50 with nl=10−2n_l=10^-2 (i.e. 1%1\% relative sensor noise). For the stochastic ADAM phase, the NtrN_tr trajectories are partitioned into NBN_B batches of size Ntr/NBN_tr/N_B, and the optimizer takes one gradient step per batch (without resampling) usually to escape poor local minima/saddle points. One pass over all the batches is an epoch. We set NB=5N_B=5 and run for 250 epochs. We freeze the batch ordering across all 20 initializations to remove shuffling as a confounder. Then, we switch to L-BFGS for 250 more epochs. Note no batch partitioning is needed for L-BFGS, since it is a deterministic optimizer. The architecture of the neural networks is also kept fixed across all cases. Unlike eqs. 45 and 46, which scale residuals by σy _y, we report a mean-normalized RMSE on the measurement and validation sets to read error as a fraction of the mean concentration, RMSEs=(1NsNt∑i=1Ns∑k=0Nt|yi(k)−y^i(k)μy|2)1/2RMSE_s= ( 1N_sN_t _i=1^N_s _k=0^N_t | y_i(k)- y_i(k) _y |^2 )^1/2 (66) for s∈train,vals∈\train,val\ where μy _y is the mean output value, and NsN_s the number of trajectories in set s; e.g. 0.020.02 means a 2%2\% error relative to the mean concentrations. (b) Black-box model We train the black-box model with the same strategy as the structured models, taking care to remove any confounding factors across the 20 initializations. We label eq. 47 as bbox and add a noisy-measurement case bbox_noise, analogous to cmeas_noise. A summary of all six cases considered in the paper is provided in table 1. We set NB=5N_B=5 and run for 250 epochs of ADAM followed by 250 epochs of L-BFGS. We use the same mean-normalized RMSE of eq. 66 to report the training and validation errors. (c) Validation The training loss, validation loss and training time for all the identified models are reported in table 3. For the structured models we first look at the quantity a practitioner can actually observe, i.e., the concentration predictions on the validation set. Figure 4 overlays these predictions on the validation measurements for all four cases. The fits are excellent, and the ensemble bands are so tight that they are barely visible. Every random initialization produces essentially the same predictions, so the 5–95% quantile range collapses onto the mean. Figure 4: Concentration predictions on the validation set for (a) rmeas, (b) cmeas_rinit, (c) cmeas and (d) cmeas_noise. Measurements are shown in red. The shaded bands are the ensemble spread, taken as the 5–95% quantiles of the 20 random initializations. In our simulation setting, we can also inspect a quantity the practitioner never sees, the plant reaction rates. The parity plots of the model rates against the truth in fig. 5 tell a very different story. The rate predictions degrade steadily from (a) rmeas to (d) cmeas_noise, even though the concentration predictions above were indistinguishable. Figure 5: Parity plots of the reaction rates [r1,r2]T[r_1,r_2]^T on the validation set for (a) rmeas, (b) cmeas_rinit, (c) cmeas and (d) cmeas_noise. Top row shows r1r_1 and bottom row r2r_2. The dotted line is ideal parity. The bands are the ensemble spread, taken as the 5–95% quantiles across the 20 random initializations. Despite discarding all physical structure, the black-box models attain a satisfactory fit when overlaid on the validation measurements, as shown in fig. 6, though less accurate than that of the structured models. Here too the ensemble bands stay tight. The black-box models thus pass every check a practitioner would run. Figure 6: Concentration predictions on the validation set for the black-box model, (a) bbox and (b) bbox_noise. Measurements are shown in red. The bands are the ensemble spread, taken as the 5–95% quantiles of the 20 random initializations. Table 3: Structured- and black-box-model results across cases. Values are reported as mean ± 1.96σ1.96\,σ over 20 random initializations (95% interval). rmsemeas and rmseval are the mean-normalized errors of eq. 66 on the measurement and validation sets; training time is wall-clock. Model Case N architecture rmsemeas rmseval Train time (hr) Structured rmeas rNN,r1=[4,64,1]r_N,r_1=[4,64,1] rNN,r2=[3,64,1]r_N,r_2=[3,64,1] 0.00067±0.000120.00067± 0.00012 0.0069±0.00210.0069± 0.0021 0.22±0.080.22± 0.08 cmeas_rinit 0.00435±0.000600.00435± 0.00060 0.0146±0.00240.0146± 0.0024 11.97±0.9411.97± 0.94 cmeas 0.00966±0.001310.00966± 0.00131 0.0219±0.00350.0219± 0.0035 12.27±1.1712.27± 1.17 cmeas_noise 0.01472±0.000640.01472± 0.00064 0.0217±0.00300.0217± 0.0030 12.45±0.9912.45± 0.99 Black-box bbox fNN=[16,128,7]f_N=[16,128,7] 0.03224±0.001870.03224± 0.00187 0.0484±0.00510.0484± 0.0051 12.01±1.0612.01± 1.06 bbox_noise 0.03410±0.001770.03410± 0.00177 0.0481±0.00570.0481± 0.0057 12.04±0.9012.04± 0.90 Herein lies the trap. Every case we considered so far reports a low training and validation RMSE and a tight ensemble spread across initializations, so nothing the practitioner can compute from the data flags a problem. Yet, as we show next, five of these six “validated” models are useless for RTO. 4.3 Optimization landscape Instead of solving the RTO problem eq. 43 over all five degrees of freedom (i.e. [NEe1,NO1,NA1,α,T]T[N_Ee^1,N_O^1,N_A^1,α,T]^T), we first solve a simpler problem over only two degrees of freedom. To see how each model reshapes the optimization landscape, we vary only the reactor inlet temperature T and the gas recycle fraction α over their box constraints [u¯,u¯][ u, u] from eq. 43, keeping the other three degrees of freedom fixed at the plant optimum. The plant profit ℓ(T,α) (T,α) is then a 2D surface we can draw. For model-based RTO, each trained model operates on its own predicted profit ℓ^(T,α) (T,α). It is distinct from ℓu^s _ u_s, the plant profit evaluated at the operating point the model recommends. This second quantity is what enters the profit-loss metric Δℓ used later. Figure 7 shows the plant landscape ℓ , and fig. 8 the model landscape ℓ that each trained model presents to the RTO optimizer, for one of the worst-performing members of the ensemble of 20, on the same axes and profit scale. Figure 7: Plant profit landscape ℓ(T,α) (T,α) with the remaining degrees of freedom fixed at the plant optimum. The star marks the plant optimum (T⋆,α⋆)(T ,α ). The surface is unimodal and well-conditioned. The plant landscape in fig. 7 has a single peak at (T⋆,α⋆)(T ,α ), consistent with the well-conditioned optimum discussed in section 2.2. Figure 8: Model-predicted profit landscape ℓ^(T,α) (T,α) that each trained model presents to the RTO optimizer. Panels, in order of degrading structure and data, are (a) rmeas, (b) cmeas_rinit, (c) cmeas, (d) cmeas_noise, (e) bbox, and (f) bbox_noise. The star marks the plant optimum as a fixed reference and the circles each model’s own optimum. Each panel shows one of the worst-performing members of the 20 initializations. Figure 8 tells a different story. Trained on rate measurements, rmeas reproduces both the plant landscape and its peak. The rmeas case also helps us diagnose that the training dataset is informative enough to identify the plant rate structure. However, cmeas_rinit (table 1) is the telling case. Initialized at the same N weights as rmeas that recovers the plant peak, the stochastic training optimizer deforms the surface so that cmeas_rinit no longer agrees with the landscape it started from. The same drift is visible in the loss curves in of the supplementary material, where ADAM moves away from the initialization and L-BFGS recovers only part of it. We trained cmeas_rinit at the same learning rate (i.e., 0.05) tuned for cmeas and cmeas_noise. A smaller learning rate does damp the drift, only because it preserves an optimum that was supplied, not one that training located. An optimizer should not move away from an optimum simply because of a user-defined hyperparameter. cmeas, cmeas_noise, bbox and bbox_noise distort the surface further with multiple local maxima. cmeas even admits an optimum in a near-saddle region, where the curvature is indefinite and gradient-correction schemes such as modifier adaptation would fail to converge. 4.4 Economic performance under RTO We next solve the actual RTO problem eq. 43 and evaluate the profit loss Δℓ over the 20 ensembles for each case of the structured and the black-box models over 3 multi-start RTO seeds. Within each multi-start seed, the RTO problem eq. 43 for each model is initialized with the same guess (xs,us)(x_s,u_s), which is obtained by perturbing the midpoint of the training data range. This ensures RTO initialization is not a confounder within a seed. As shown in fig. 9, rmeas recovers the plant optimum with negligible profit loss and no variability across the multi-start seeds. This result again confirms that the training dataset is informative enough to identify the plant rate structure. Figure 9: Profit loss Δℓ for (a) rmeas, (b) cmeas_rinit, (c) cmeas and (d) cmeas_noise. Columns are cases and rows are the 3 multi-start RTO seeds; each panel is the distribution over the 20 ensemble members, plotted as relative frequency. The dotted line marks the best achievable profit loss, 0%0\%. Every other case, however, incurs a significant profit loss that varies across the RTO seeds. The worst-case profit loss for structured models is just over 30%30\% as seen from fig. 9 panel (d) for seed 1. Models that were visually indistinguishable on the validation dataset predict multiple optima, even though the same multistart search returns a single optimum for the plant. This disparity among RTO solutions is not an artifact of noise, since apart from cmeas_noise all results for structured models used noise-free measurements. Nor would training the models longer help, since the validation loss for the models has already plateaued within the chosen epoch budget. As noted previously, a good starting point does not help either. Even after warm-starting from the rmeas weights and biases that recover the plant optimum, cmeas_rinit drifts and returns a near-zero loss only occasionally. A natural objection, often raised in support of black-box models, is that imposing structure itself is the culprit. By fixing the stoichiometry with mass balances and thermodynamics, we may have imposed a model class too rigid for the plant, and letting the data determine the entire dynamics might do better on RTO. We tested the claim directly with the black-box models. Their economic performance, however, is still worse for the actual RTO problem. The profit loss in fig. 10 exceeds that of every structured case, with the worst-case exceeding 50%50\% as seen from panel (b) for seed 3. Figure 10: Profit loss Δℓ for (a) bbox and (b) bbox_noise. Columns are cases and rows are the 3 multi-start RTO seeds; each panel is the distribution over the 20 ensemble members, plotted as relative frequency. The dotted line marks the best achievable profit loss, 0%0\%. 5 Conclusions In this paper, we tested whether a data-driven model that predicts plant measurements accurately can be trusted to drive economic decisions. Using the vinyl acetate process as a benchmark, we trained hybrid (structured) models, which embed a neural-network closure for reaction rates in the known mass balances and thermodynamics, and a fully data-driven black-box model. Within each model class, we considered several cases that differ in the measurements used for training and the noise level. We tested each case in real-time optimization. Across the ensemble based on random initializations, each model and its variants cleared the usual validation checks, with accurate fits and tight ensemble spread on the validation dataset, yet their RTO solutions showed a wide disparity. In the worst case, the structured models lost 30%30\% of profit and the black-box cases over 50%50\%. With losses this large, we may do better for economic profitability by not using real-time optimization at all. Additionally, under the same multistart search that returns a single optimum for the plant, the models returned many spurious RTO solutions. We ran these simulations under the most ideal settings we could construct. We took care to have the measurements at the plant optimum bracketed within the training dataset, ensuring that the models are tested purely on interpolation. We trained the models on both noise-free and noisy measurements. In fact, in the cmeas_rinit case the N weights and biases that recover the plant RTO landscape were handed to the training optimizer as its starting point. Training on concentration data drove those parameters away to a suboptimal fit whose RTO landscape no longer matched the one it started from. The observed optimizer drift is not a numerical failure. Stochastic optimizers such as ADAM, which are workhorses of ML/AI frameworks, advance without a line search, and that is precisely what lets them escape poor local minima. The same property, however, can let them drift away from a good basin. Switching to a deterministic optimizer, L-BFGS, for the later epochs does not repair the drift; L-BFGS performs a line search but by the time it takes over, the initial ADAM phase has already left the good basin it was initialized in, so L-BFGS only refines within whichever basin it inherits. A practitioner does not expect an optimizer to leave an optimal solution it was handed. The effect does not show up in the training or validation fit, but in the RTO decision. Nor is starting L-BFGS from the correct basin a remedy, since it presumes one already holds the RTO-optimal parameters. This fragility of the standard ML training recipe is what we highlight. This fragility can only worsen in the presence of process noise, unmodeled states or disturbances, or insufficiently exciting closed-loop data. We ran these investigations with the best software and tuned hyperparameters, and the models still failed under these ideal conditions. Any industrial vendor proposing to solve a similar problem should first demonstrate recovery of the optimum on a benchmark of the type proposed here. Short of that, the attempt remains unreliable. Additional regularization of the structured models could help overcome such fragility. Future work could test whether imposing monotonicity constraints on certain inputs (e.g., temperature) in the N-approximated rate laws, to match insights from process chemistry, improves RTO performance. Another avenue is to impose a Lipschitz bound on the neural-network curvature and test whether this makes data-driven model-based RTO more robust. Acknowledgment The funding for this work was provided by the Dow Chemical Company. PD acknowledges support from the Mitsubishi Chemical Graduate Fellowship. Use was made of computational facilities purchased with funds from the National Science Foundation (CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF DMR 2308708) at UC Santa Barbara. Data availability The code and simulated data used to reproduce all results are available at https://dakeprithvi.github.io/2026c_structure_id_plant and Dake et al. (2026). Appendix A Process units Every separator and splitter acts as a linear split: an inlet molar flow NjinN_j in of species j divides into a split-off stream and a remainder, Njsplit=δjNjin,Njrem=(1−δj)Njin,j∈J,N_j split= _j\,N_j in, 18.49988ptN_j rem=(1- _j)\,N_j in, 18.49988ptj∈ J, where J=Ee,Ea,A,V,W,O,CJ=\Ee,Ea,A,V,W,O,C\ denotes ethylene, ethane, acetic acid, vinyl acetate, water, oxygen, and carbon dioxide. Table 4 lists each unit, its inlet/outlet streams, and its split fraction. Each unit carries its own symbol for δj _j so the flowsheet equations stay unambiguous, ϕj _j for the flash drum, γj _j for the distillation column, μj _j for the absorber, and the recycle fractions α and β for splitters I and I. Table 4: Linear split for each separation unit and splitter. Stream superscripts follow the flowsheet; the last column gives each unit’s split fraction, the fraction of NjinN_j in routed to the split-off stream. Unit In Split/Rem Split fraction Flash drum N3N^3 N5/N4N^5/N^4 ϕj=1 _j=1: Ee,Ea,O,C; 00: A,V,W Distillation column N4N^4 N8/N10N^8/N^10 γj=1 _j=1: A; 00: otherwise Absorber N5N^5 N11/N6N^11/N^6 μj=0.03 _j=0.03: C; 00: otherwise Splitter I N6N^6 N7/N9N^7/N^9 α (all j) Splitter I N8N^8 N13/N12N^13/N^12 β (all j) Appendix B Parameter values The price-vector is given as, (p1,p2,p3,p4,Ee,p4,O,p4,A,p5,p6)=(p_1,p_2,p_3,p_4,Ee,p_4,O,p_4,A,p_5,p_6)= (300, 33, 50, 100, 50, 50, 30, 0) $/mol \$mol Table 5: Parameter values used for the VAc model. Parameter Value Units k1k_1 1.312×10−21.312× 10^-2 mol/kgs(m3/mol)3.2 molkg\,s ( m^3mol )^3.2 E1E_1 40500 J/mol k2k_2 2.3632.363 mol/kgs(m3/mol)2 molkg\,s ( m^3mol )^2 E2E_2 52500 J/mol R 8.314 J/molK JmolK P 882.529 kPa VRV_R 6.68 m3m^3 ϵε 0.8 - ρc _c 385 kg/m3 kgm^3 λ 0.001 - u¯ u (10.88, 7.8, 12.01, 0.4, 0.6, 384.39) (mol/s mols, mol/s mols, mol/s mols, -, -, KK) u¯ u (20.20, 14.49, 22.30, 0.74, 0.6, 469.80) (mol/s mols, mol/s mols, mol/s mols, -, -, KK) sus_u (13.85, 13.85, 13.08, 1, 1, 423.15) (mol/s mols, mol/s mols, mol/s mols, -, -, KK) References Andersson et al. (2019) J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl CasADi—a software framework for nonlinear optimization and optimal control. Math. Prog. Comp. 11 (1), p. 1–36. External Links: Document Cited by: §2.2. Bindlish (2025) R. Bindlish An Industrial Perspective on Applying Process Knowledge and Mathematical Programming for RTO, MPC, NMPC, and Power Scheduling. Ind. Eng. Chem. Res. 64 (50), p. 23810–23823. Cited by: §1.1, §2.3. Bonvin and Pannocchia (2024) D. Bonvin and G. Pannocchia On speeding-up modifier-adaptation schemes for real-time optimization. Comput. Chem. Eng. 191, p. 108839. Cited by: §1. Box (1957) G. E. P. Box Evolutionary operation: a method for increasing industrial productivity. Journal of the Royal Statistical Society Series C 6 (2), p. 81–101. External Links: Document, Link Cited by: §1. Bradbury et al. (2018) J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang JAX: composable transformations of Python+NumPy programs. External Links: Link Cited by: §3.2. Câmara et al. (2016) M. M. Câmara, A. D. Quelhas, and J. C. Pinto Performance evaluation of real industrial RTO systems. Processes 4 (4), p. 44. Cited by: §1. Chachuat et al. (2009) B. Chachuat, B. Srinivasan, and D. Bonvin Adaptation strategies for real-time optimization. Comput. Chem. Eng. 33 (10), p. 1557 – 1567. Cited by: §1. Chen et al. (2018) R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud Neural ordinary differential equations. Adv. Neural Inf. Process. Syst. 31, p. 6571–6583. Cited by: §1.1. Chen et al. (2003) R. Chen, K. Dave, T. J. McAvoy, and M. Luyben A nonlinear dynamic model of a vinyl acetate process. Ind. Eng. Chem. Res. 42 (20), p. 4478–4487. Cited by: §2. Constante et al. (2026) G. E. Constante, H. Chen, and C. Li Enforcing hard linear constraints in deep learning models with decision rules. Vol. 38, p. 105075–105100. Cited by: §1.1. Crowe (1996) C. M. Crowe Data reconciliation—progress and challenges. J. Proc. Cont. 6 (2-3), p. 89–98. Cited by: §1. Cutler and Perry (1983) C. R. Cutler and R. T. Perry Real time optimization with multivariable control is required to maximize profits. Comput. Chem. Eng. 7, p. 663–667. Cited by: §1. Dake et al. (2026) P. Dake, R. Bindlish, and J. Rawlings Dataset for ‘A tale of perfect fit and phantom optima: how data-driven models can fail in real- time optimization’. Zenodo. External Links: Document, Link Cited by: Data availability. Dake et al. (2024) P. Dake, M. R. Ilagan, S. Banerjee, S. L. Scott, and J. B. Rawlings Identifying kinetic models from reactor measurements. J. Vac. Sci. Technol. A. 42 (6). External Links: Document, ISSN 1520-8559 Cited by: §4.1. Darby et al. (2011) M. L. Darby, M. Nikolaou, J. Jones, and D. Nicholson RTO: An overview and assessment of current practice. J. Proc. Cont. 21 (6), p. 874–884. External Links: ISSN 0959-1524 Cited by: §1, §1. Forbes et al. (1994) J. F. Forbes, T. E. Marlin, and J. F. MacGregor Model adequacy requirements for optimizing plant operations. Comput. Chem. Eng. 18 (6), p. 497–510. External Links: ISSN 0098-1354 Cited by: §1.1. Fronk and Petzold (2025) C. Fronk and L. Petzold The vanishing gradient problem for stiff neural differential equations. Chaos 35 (11). Cited by: §3.3. Golder et al. (2025) R. Golder, B. N. Roy, and M. Hasan DAE-hardnet: a physics constrained neural network enforcing differential-algebraic hard constraints. Cited by: §1.1. Heek et al. (2024) Flax: a neural network library and ecosystem for JAX External Links: Link Cited by: §3.2. Hornik et al. (1989) K. Hornik, M. Stinchcombe, and H. White Multilayer feedforward networks are universal approximators. Neural Netw. 2 (5), p. 359–366. Cited by: §1.1. Kidger (2021) P. Kidger On Neural Differential Equations. Ph.D. Thesis, University of Oxford. Cited by: §3.2. Kingma and Ba (2014) D. P. Kingma and J. Ba Adam: a method for stochastic optimization. arXiv preprint arXiv:1412.6980. Cited by: §3.2. Krishnamoorthy et al. (2018) D. Krishnamoorthy, B. Foss, and S. Skogestad Steady-state real-time optimization using transient measurements. Comput. Chem. Eng. 115, p. 34–45. Cited by: §1. Krishnamoorthy and Skogestad (2022) D. Krishnamoorthy and S. Skogestad Real-time optimization as a feedback control problem–a review. Comput. Chem. Eng. 161, p. 107723. Cited by: §1. Kumar and Rawlings (2023) P. Kumar and J. B. Rawlings Structured nonlinear process modeling using neural networks and application to economic optimization. Comput. Chem. Eng. 177, p. 108314. External Links: Document Cited by: §1.1, §2.3, §3.1. LeCun et al. (2002) Y. LeCun, L. Bottou, G. B. Orr, and K. Müller Efficient backprop. In Neural networks: Tricks of the trade, p. 9–50. Cited by: Table 1. Liu and Nocedal (1989) D. C. Liu and J. Nocedal On the limited memory bfgs method for large scale optimization. Math. Prog. 45 (1), p. 503–528. Cited by: §3.2. Lueg et al. (2025) L. R. Lueg, V. Alves, D. Schicksnus, J. R. Kitchin, C. D. Laird, and L. T. Biegler A simultaneous approach for training neural differential-algebraic systems of equations. arXiv preprint arXiv:2504.04665. Cited by: §1.1. Luyben and Tyréus (1998) M. L. Luyben and B. D. Tyréus An industrial design/control study for the vinyl acetate monomer process. Comput. Chem. Eng. 22 (7-8), p. 867–877. Cited by: §2. Mukherjee and Zavala (2026) A. Mukherjee and V. M. Zavala Physics-constrained machine learning for chemical engineering. Curr. Opin. Chem. Eng. 51, p. 101228. Cited by: §1.1, §1.1. Naysmith and Douglas (1995) M. R. Naysmith and P. L. Douglas Review of real time optimization in the chemical process industries. 3 (2), p. 67–87. Cited by: §1. Pantelides et al. (2025) C. Pantelides, M. Baldea, A. T. Georgiou, B. Gopaluni, M. Mehmet, K. Sheth, V. M. Zavala, and C. Georgakis From Automated to Autonomous Process Operations. Comput. Chem. Eng. 196. External Links: Document Cited by: §1.1. Papasavvas et al. (2019) A. Papasavvas, T. de Avila Ferreira, A. G. Marchetti, and D. Bonvin Analysis of output modifier adaptation for real-time optimization. Comput. Chem. Eng. 121, p. 285–293. Cited by: §1. Peng et al. (2025) W. Peng, Z. Liu, and M. Boy Spin-ode: stiff physics-informed neural ode for chemical reaction rate estimation. Cited by: §3.3. Raissi et al. (2019) M. Raissi, P. Perdikaris, and G. E. Karniadakis Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 378, p. 686–707. Cited by: §1.1. Rawlings et al. (2012) J. B. Rawlings, D. Angeli, and C. Bates Fundamentals of economic model predictive control. In IEEE Conference on Decision and Control (CDC), Maui, HI, p. 3851–3861. Cited by: §1. Roberts and Williams (1981) P. D. Roberts and T. W. C. Williams On an algorithm for combined system optimisation and parameter estimation. Automatica 17 (1), p. 199–209. Cited by: §1. Sansana et al. (2021) J. Sansana, M. N. Joswiak, I. Castillo, Z. Wang, R. Rendall, L. H. Chiang, and M. S. Reis Recent trends on hybrid modeling for industry 4.0. Comput. Chem. Eng. 151, p. 107365. Cited by: §1.1, §1.1. Schweidtmann et al. (2024) A. M. Schweidtmann, D. Zhang, and M. Von Stosch A review and perspective on hybrid modeling methodologies. Digital Chemical Engineering 10, p. 100136. Cited by: §1.1, §1.1. Seborg et al. (2017) D. E. Seborg, T. F. Edgar, D. A. Mellichamp, and F. J. Doyle Process dynamics and control. fourth edition, John Wiley and Sons, New York. Cited by: §1. Sholokhov et al. (2023) A. Sholokhov, Y. Liu, H. Mansour, and S. Nabi Physics-informed neural ode (pinode): embedding physics into models using collocation points. Sci. Rep. 13 (1), p. 10166. Cited by: §1.1. Skogestad (2000) S. Skogestad Self-optimizing control: the missing link between steady-state optimization and control. Comput. Chem. Eng. 24 (2-7), p. 569–575. Cited by: §1. Thompson et al. (2026) J. Thompson, B. M. Connors, V. M. Zavala, and O. S. Venturelli Physics-constrained neural ordinary differential equation models to discover and predict microbial community dynamics. Proc. Natl. Acad. Sci. USA 123 (13), p. e2517661123. Cited by: §1.1. Wächter and Biegler (2006) A. Wächter and L. T. Biegler On the implementation of a primal-dual interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Prog. 106 (1), p. 25–57. Cited by: §2.2. Ward et al. (2004) J. D. Ward, D. A. Mellichamp, and M. F. Doherty Importance of process chemistry in selecting the operating policy for plants with recycle. Ind. Eng. Chem. Res. 43 (14), p. 3957–3971. Cited by: §2. Supplementary material for “A tale of perfect fit and phantom optima: how data-driven models can fail in real-time optimization” S1 Learning-rate sweep and optimizer comparison For reference, we restate the training problems whose learning rate is swept here. The concentration-measurement cases cmeas_rinit, cmeas, and cmeas_noise solve the structured training problem, Eq. (9) of the main paper, ΘN⋆ _N =argminΘN(1NtrNt∑i=1Ntr∑k=0Nt|yi(k)−y^i(k)σy|2)1/2 = _ _N ( 1N_trN_t _i=1^N_tr _k=0^N_t | y_i(k)- y_i(k) _y |^2 )^1/2 s.t.x^˙i s.t. x_i =f(x^i,ui,ΘN) =f( x_i,u_i; _N) y^i y_i =h(x^i,ui) =h( x_i,u_i) x^i(0) x_i(0) =f0(yi(0)) =f_0(y_i(0)) (1) while the two black-box cases bbox and bbox_noise solve the black-box training problem, Eq. (11) of the main paper, θfNN⋆ θ _f_N =argminθfNN(1NtrNt∑i=1Ntr∑k=0Nt|yi(k)−y^i(k)σy|2)1/2 = _ _f_N ( 1N_trN_t _i=1^N_tr _k=0^N_t | y_i(k)- y_i(k) _y |^2 )^1/2 s.t.x^˙i s.t. x_i =fNN(x^i,ui,θfNN) =f_N( x_i,u_i; _f_N) y^i y_i =h(x^i,ui) =h( x_i,u_i) x^i(0) x_i(0) =f0(yi(0)) =f_0(y_i(0)) (2) The symbols are defined in the main paper. Every sweep in this section is performed at the fixed batch count NB=5N_B=5 used throughout the paper (the NtrN_tr training trajectories are partitioned into NBN_B batches of size Ntr/NBN_tr/N_B during the stochastic ADAM phase); only the learning rate is varied. Figures S1–S5 show the learning-rate sweep for each case, and every figure has the same four-panel layout. The top row is the training loss and the bottom row the validation loss, both on a logarithmic scale; the left column is ADAM alone and the right column is ADAM followed by L-BFGS. Panel (a) is therefore the training loss under ADAM, (b) the training loss under ADAM and then L-BFGS, (c) the validation loss under ADAM, and (d) the validation loss under ADAM and then L-BFGS. Within each panel we overlay one curve per learning rate across the full epoch budget, so a single panel shows both the learning-rate sweep and whether the loss has plateaued within that budget. In the right column, panels (b) and (d), the dashed vertical line at epoch 250 marks the switch from Adam to L-BFGS. For cmeas and cmeas_noise, and for the two black-box cases, we select the learning rate that gives the lowest validation loss in panel (d) and confirm from the same panel that the validation loss has plateaued within the epoch budget. cmeas_rinit is read differently. It is initialized at the weights at which the model recovers the plant optimum, so no learning rate is selected from it. In the paper, however we show results for the same learning rate as chosen for cmeas and cmeas_noise, i.e., 0.05. The dash-dotted horizontal line, labelled learning rate 00, is the loss at that initialization before any training step is taken. Figure S1: cmeas_rinit: drift away from the initialization. The model is initialized at the weights at which it recovers the plant optimum, marked by the dash-dotted line (learning rate 00, no training step taken). Under ADAM (left) the loss moves away from that line as training proceeds, so the optimizer leaves the basin it was handed. After the switch to L-BFGS at epoch 250 (right) the loss returns toward the initialization, so L-BFGS recovers part of what ADAM gave up. As the learning rate decreases, ADAM takes smaller steps and the curves stay at the initialization; this preserves the optimum only because it was supplied, not because training located it. Figure S2: Learning-rate sweep for cmeas. The chosen learning rate is 0.050.05. Figure S3: Learning-rate sweep for cmeas_noise. The chosen learning rate is 0.050.05. Figure S4: Learning-rate sweep for bbox. The chosen learning rate is 0.010.01. Figure S5: Learning-rate sweep for bbox_noise. The chosen learning rate is 0.010.01.