Catalyst-Activity-Aware Self-Optimizing Digital Twin for o-Xylene Oxidation to Phthalic Anhydride in a Catalyst-Deactivating Fixed-Bed Reactor

1. Introduction

Catalytic fixed-bed reactors are central to many chemical and petrochemical processes, including oxidation, hydrogenation, reforming, desulfurization, and synthesis reactions. Their performance depends strongly on catalyst activity, temperature distribution, feed composition, heat-transfer behavior, and operating constraints. Catalyst deactivation is a persistent limitation in heterogeneous catalytic processes because activity and selectivity can decline through poisoning, fouling, coking, sintering, phase transformation, or loss of active surface area [1,2,3]. In industrial operation, reactor targets are often selected from design calculations, commissioning studies, offline optimization, or periodic plant tests. These targets may remain fixed for long periods, even though the reactor itself changes continuously as the catalyst deactivates and process conditions drift.

Catalyst deactivation is particularly important in highly exothermic fixed-bed reactors. A reduction in catalyst activity may lower conversion, alter selectivity, shift hot-spot locations, and change the relationship between manipulated variables and product yield. Operators may compensate by increasing inlet temperature, changing coolant conditions, or modifying feed rates [4,5,6]. Without a reliable estimate of the current catalyst condition, however, such adjustments may be conservative, delayed, or based on indirect operating experience rather than systematic optimization.

Digital twins provide a promising route for improving reactor operation because they combine physics-based process models, real-time measurements, and computational decision-making [7,8,9,10]. In the process industries, digital-twin implementations are increasingly discussed for monitoring, prediction, predictive maintenance, process optimization, and operational decision support, but deployment remains challenging because models must remain synchronized with changing plant behavior and uncertain measurements [11]. Many existing digital twins are used primarily for monitoring, fault detection, or operator support. A more powerful role emerges when the digital twin can estimate hidden process states and propagate those estimates into operating decisions. For catalytic reactors, one of the most consequential hidden states is the catalyst activity profile. Since catalyst activity is not directly measured during normal operation, it must be inferred from available process measurements such as axial temperature, outlet composition, flow rate, and pressure.

This work develops and evaluates, through simulation, a state-estimation-based self-optimizing digital twin framework for catalyst-deactivating fixed-bed reactors. The key idea is to treat catalyst activity as an estimated hidden state rather than as a fixed model parameter. Process measurements and prior information are reconciled through a constrained moving-window estimation problem with physical bounds on catalyst activity [12,13,14]. The estimated catalyst condition is then passed to an optimization layer that updates reactor operating targets. This structure connects hidden-state estimation with self-optimizing operation, in which operating targets are selected to maintain near-optimal performance despite disturbances and process changes [15,16,17]. The resulting digital twin is therefore evaluated as an adaptive decision-support framework rather than solely as a prediction tool.

The proposed framework is demonstrated using the oxidation of o-xylene to phthalic anhydride in a heat-exchanged fixed-bed reactor. This system is suitable because catalyst activity affects both the axial temperature profile and outlet conversion, providing measurable information about deactivation. It is also a representative reaction-engineering problem because the desired partial oxidation competes with over-oxidation pathways under strong heat release and hot-spot constraints [5,6,18,19]. The case study combines a reactor model benchmarked against published results, axial activity-profile estimation from sparse synthetic measurements, adaptive thermal-target optimization, and robustness tests under different noise levels, sensor configurations, regularization choices, and kinetic mismatch.

Catalyst deactivation in fixed-bed reactors has been studied from several complementary perspectives. One group of studies focuses on physical or semi-empirical deactivation modeling, where loss of catalyst activity is represented through temperature-, composition-, or time-dependent deactivation laws. For o-xylene oxidation to phthalic anhydride, industrial deactivation has been analyzed using temperature-profile data and detailed reactor models to infer how the catalyst condition changes along the bed [4]. More general activity-profile modeling methods have also been proposed to approximate spatially distributed activity loss using reduced-order representations [20]. These studies are valuable because they clarify how catalyst aging affects reactor behavior, but their primary objective is deactivation description or model reconstruction rather than online operating-target adaptation.

A second group of studies addresses catalyst activity estimation. In this direction, activity profiles are inferred from process measurements such as axial temperature and outlet composition. The work of Cheng et al. is particularly relevant because it demonstrated that temperature, composition, and catalyst activity profiles can be estimated in fixed-bed reactors with decaying catalysts using process measurements and model-based estimation [21]. This establishes the feasibility of reconstructing hidden catalyst states from measurable reactor outputs. However, activity estimation is often used mainly for monitoring, diagnosis, or model updating, and the estimated activity is not necessarily propagated into an optimization layer that changes reactor targets.

A third line of work focuses on reactor optimization, optimal temperature policies, and digital or model-based reactor design. For phthalic anhydride synthesis, reactor optimization studies have investigated how temperature policy and reactor design influence selectivity, yield, and safe operation [5]. More recent model-based digital design concepts for fixed-bed catalytic reactors emphasize the role of mechanistic models in design, scale-up, optimization, and possible online implementation [22]. These optimization studies are important, but they commonly assume that catalyst activity is known, fixed, or updated outside the optimization loop.

Taken together, the literature provides strong foundations for deactivation modeling, catalyst activity estimation, and reactor optimization. However, these elements are commonly studied separately. Activity estimates are often used for monitoring, diagnosis, or model reconstruction without being propagated into operating-target optimization, whereas reactor optimization studies commonly assume that catalyst activity is known, fixed, or updated outside the optimization loop. To the authors’ knowledge, limited prior work has examined the complete pathway from sparse reactor measurements to spatial activity estimation and then to catalyst-aware target adaptation. This work addresses that integration gap through a simulation-based proof of concept linking physics-based reactor prediction, moving-window constrained activity estimation, and self-optimizing target updates.

Accordingly, the main research question is whether catalyst activity inferred sequentially from sparse process measurements can be used within a digital twin to recover reactor performance as deactivation progresses, and how sensitive the resulting operating decisions are to measurement quality, sensor configuration, estimator structure, and model–plant mismatch.

The paper makes four principal contributions:

• A general digital-twin workflow is formulated in which an estimated catalyst state is used directly in operating-target optimization.
• The workflow is instantiated for o-xylene oxidation using a literature-benchmarked heat-exchanged fixed-bed reactor model with spatially nonuniform catalyst activity.
• Sparse axial temperatures and outlet conversion are combined in a constrained moving-window estimator, and the propagation of estimation error into optimized operating targets is quantified.
• The workflow is evaluated using repeated noise realizations, sensor-density tests, a uniform-activity comparator, regularization ablation, and kinetic model–plant mismatch.

The novelty therefore lies in the demonstrated integration and evaluation of the closed decision pathway, rather than in proposing a new kinetic model or claiming the first use of catalyst activity estimation.

Table 1 positions the present work relative to representative studies across six relevant themes. The contribution is not the reactor model or activity estimator in isolation, but their integration so that an estimated catalyst condition becomes a direct input to an operating-target update layer.

The remainder of the paper is organized as follows. Section 2 presents the proposed methodology, including the closed-loop workflow, activity estimation problem, self-optimizing operating layer, and robustness evaluation procedure. Section 3 describes the o-xylene oxidation case study, reactor model, kinetic expressions, assumptions, validation basis, and catalyst activity representation. Section 4 presents and discusses the results, including activity observability, estimation performance, self-optimizing operation, robustness tests, novelty positioning, and limitations. Section 5 summarizes the conclusions.

2. Methodology

2.1. General Self-Optimizing Digital Twin Framework

The proposed framework provides a general structure for catalyst-deactivating fixed-bed reactors. The central idea is to combine process measurements, a physics-based reactor digital twin, a sequential catalyst activity estimator, and a self-optimizing operating layer in a catalyst-state-aware decision workflow. Figure 1 illustrates the intended closed-loop architecture: process measurements are compared with model predictions to estimate the hidden catalyst activity profile, and the estimate is then used to update reactor operating targets. In the present proof-of-concept implementation, estimation is first evaluated along a prescribed deactivation trajectory at fixed nominal inputs, after which the estimated activity sequence is supplied to the target-optimization layer. Thus, the study demonstrates the complete measurement–estimation–optimization decision pathway, but it does not simulate actuator dynamics or feed each optimized target back into the synthetic plant measurement generator at the next sampling instant.

The framework is applicable to catalytic reactor configurations where catalyst degradation affects measurable outputs such as temperature, conversion, selectivity, or yield. Representative examples include multi-bed adiabatic reactors, heat-exchanged fixed-bed reactors, and staged catalytic systems. The framework consists of four interacting layers:

1.

Physical reactor and measurement system: provides process measurements such as axial temperature, outlet composition, flow rate, and pressure.

2.

Physics-based digital twin: predicts reactor states from a mechanistic model given current inputs and the estimated catalyst activity.

3.

Sequential catalyst activity estimator: uses the mismatch between measured and predicted outputs to infer the hidden catalyst activity profile at each sampling step.

4.

Self-optimizing operating layer: uses the estimated activity to compute updated reactor operating targets subject to process constraints.

The general closed-loop workflow is

$$\begin{array}{c}\mathrm{reactor}\to \mathrm{measurements}\to \mathrm{digital}\mathrm{twin}\hfill \\ \hfill \to \mathrm{state}\mathrm{estimation}\to \mathrm{self}-\mathrm{optimizing}\mathrm{layer}\to \mathrm{updated}\mathrm{targets}.\end{array}$$

(1)

This structure enables adaptive operation because the optimization problem is solved using the current estimated catalyst condition rather than assuming fixed or fresh catalyst activity.

2.1.1. Physics-Based Reactor Model

For a general catalytic reaction network containing $N_r$ reactions, reaction q may be written as

$${\displaystyle \sum _{j=1}^{N_c}}{\nu }_{jq}^{-}{A}_j\to {\displaystyle \sum _{j=1}^{N_c}}{\nu }_{jq}^{+}{A}_j,$$

(2)

where $A_j$ denotes chemical species j, $N_c$ is the number of components, and ${\nu }_{jq}^{-}$ and ${\nu }_{jq}^{+}$ are the reactant and product stoichiometric coefficients. Defining the signed coefficient ${\nu }_{jq}={\nu }_{jq}^{+}-{\nu }_{jq}^{-}$, the balance for component j in catalyst bed or reactor segment i may be expressed in catalyst-weight coordinates as

$${\displaystyle \frac{d{F}_{j,i}}{d{W}_i}}={\displaystyle \sum _{q=1}^{N_r}}{\nu }_{jq}{r}_{q,i}.$$

(3)

For an adiabatic catalyst bed, the energy balance can be written as

$${\displaystyle \frac{d{T}_i}{d{W}_i}}={\displaystyle \frac{-{\sum }_{q=1}^{N_r}\Delta H_q\left(T_i\right)r_{q,i}}{{\sum }_{j=1}^{N_c}{F}_{j,i}{C}_{p,j}\left(T_i\right)}}.$$

(4)

For heat-exchanged reactors, the energy balance is modified by adding a heat-removal term that depends on the heat-transfer coefficient, heat-transfer area, and coolant temperature.

For a key reactant A, the bed conversion is

$$X_i={\displaystyle \frac{F_{A,in,i}-F_{A,out,i}}{F_{A,in,i}}},$$

(5)

and the overall reactor conversion is

$$X_{total}={\displaystyle \frac{F_{A,0}-F_{A,out}}{F_{A,0}}}.$$

(6)

For reversible systems, equilibrium relations may provide additional feasibility constraints. The present o-xylene oxidation case is modeled using the irreversible reaction network described in Section 3.3, so no equilibrium-conversion constraint is imposed.

2.1.2. Catalyst Activity-Modified Kinetics

Catalyst activity is introduced as a slowly varying state that modifies the apparent reaction rates. For reaction q in reactor segment i, the rate is written as

$$r_{q,i}\left(t\right)=a_i\left(t\right)r_{q,i}^{int}(T_i,{\mathit{C}}_i),$$

(7)

where $r_{q,i}^{int}$ is the intrinsic rate under reference (fresh-catalyst) conditions, ${\mathit{C}}_i$ is the concentration vector, and $a_i\left(t\right)\in (0,1]$ is the local catalyst activity factor. In general, the activity factor satisfies

$$0<a_i\left(t\right)\le 1,$$

(8)

where $a_i=1$ corresponds to fresh catalyst and $a_i\to 0$ represents complete local deactivation. In practice a small positive lower bound is imposed to maintain numerical conditioning of the reactor model; the value used in the present case study is given in Section 2.3. For spatially distributed catalyst beds, activity may be represented as a continuous axial profile $a(z,t)$ or approximated by a finite number of activity nodes:

$$\mathit{a}\left(t\right)={[a_1\left(t\right),a_2\left(t\right),\dots ,a_m\left(t\right)]}^T.$$

(9)

2.1.3. Augmented State Estimation

A general process state vector may be defined as

$$\mathit{x}={[\mathit{T},\mathit{F}]}^T,$$

(10)

where $\mathit{T}$ collects the axial temperature profile and $\mathit{F}$ the species molar flow rates; derived quantities such as conversion, selectivity, and reaction rates are outputs of the model rather than independent states. To include catalyst activity, the augmented state vector is

$${\mathit{x}}_a={[\mathit{x},\mathit{a}]}^T.$$

(11)

For a dynamic implementation, the digital twin may predict process evolution using

$${\widehat{\mathit{x}}}_{k+1}=f({\widehat{\mathit{x}}}_k,{\mathit{u}}_k,{\widehat{\mathit{a}}}_k),$$

(12)

and the measurement model is

$${\widehat{\mathit{y}}}_k=h({\widehat{\mathit{x}}}_k,{\mathit{u}}_k,{\widehat{\mathit{a}}}_k).$$

(13)

For reactors where deactivation is slow relative to the process time constants, the state propagation simplifies to a quasi-steady solve ${\widehat{\mathit{x}}}_k=f({\mathit{u}}_k,{\widehat{\mathit{a}}}_k)$ at each estimation instant. This simplification is applied in the present case study and justified in Section 3.5.

2.2. Closed-Loop Workflow Algorithm

Algorithm 1 summarizes the framework-level computational workflow. The same logic can be applied to other catalyst-deactivating reactors when the activity state is observable from available measurements. In the present PA case study, the measured outputs are sparse axial temperatures and outlet conversion, the estimated hidden state is the axial catalyst activity profile, and the optimized targets are the gas and coolant inlet temperatures. As noted above, the current implementation evaluates the estimation sequence at nominal inputs and subsequently evaluates the target updates; Step 6 represents the feedback action required in a future plant-connected or fully coupled closed-loop implementation.

Algorithm 1 State-estimation-based self-optimizing digital twin workflow

Require: Physics-based reactor model $f\left(·\right)$ and measurement model $h\left(·\right)$ Require: Initial prior activity estimate ${\widehat{\mathit{a}}}_{-1}$, input bounds ${\mathit{u}}_{min},{\mathit{u}}_{max}$, and safety limit $T_{max}^{safe}$ Require: Measurement covariance $\mathit{R}$, regularization weights ${\lambda }_s,{\lambda }_p$, and sampling index $k=0,1,\dots ,N-1$ 1: for each sampling instant k do 2:     Collect process measurements ${\mathit{y}}_k$, including sparse reactor temperatures and outlet conversion. 3:     Simulate the digital twin using current inputs ${\mathit{u}}_k$ and a trial activity profile $\widehat{\mathit{a}}$. 4:     Estimate the catalyst activity profile by solving $${\widehat{\mathit{a}}}_k=arg\underset{0.05\le \widehat{\mathit{a}}\le 1}{min}{∥{\mathit{y}}_k-h({\widehat{\mathit{x}}}_k,{\mathit{u}}_k,\widehat{\mathit{a}})∥}_{{\mathit{R}}^{-1}}^2+{\lambda }_s{∥\mathit{D}\widehat{\mathit{a}}∥}_2^2+{\lambda }_p{∥\widehat{\mathit{a}}-{\widehat{\mathit{a}}}_{k-1}∥}_2^2.$$ 5:     Solve the catalyst-aware target optimization problem $${\mathit{u}}_{k+1}^{*}=arg\underset{{\mathit{u}}_{min}\le \mathit{u}\le {\mathit{u}}_{max}}{max}J({\widehat{\mathit{x}}}_k,{\widehat{\mathit{a}}}_k,\mathit{u}),$$ where J maximizes PA yield with move penalties, and any candidate violating $T_{max}>T_{max}^{safe}$ is excluded (Section 2.4). 6:     Apply or recommend the updated operating targets ${\mathit{u}}_{k+1}^{*}$; in the present study, record them for performance evaluation. 7:     Store ${\widehat{\mathit{a}}}_k$, ${\mathit{u}}_{k+1}^{*}$, predicted outputs, and performance indicators for monitoring. 8: end for

In the present implementation, $\mathit{R}=\mathrm{diag}({\sigma }_T^2{\mathit{I}}_s,{\sigma }_X^2)$ is the diagonal measurement-noise covariance, where s is the number of temperature sensors. The detailed residual construction and regularization terms are given in Section 2.3.2.

2.3. Moving-Window Constrained Activity Estimation

2.3.1. Measurement Vector

The estimator uses a measurement vector consisting of sparse axial gas temperatures and outlet conversion:

$${\mathit{y}}_k={[T(z_1,t_k),T(z_2,t_k),\dots ,T(z_s,t_k),X_{out}\left(t_k\right)]}^T,$$

(14)

where s is the number of temperature sensors. In the base case, temperature measurements are taken at four axial locations: $z=[0.5,1.0,1.5,2.0]$ m. For the robustness scenarios, the denser six-sensor layout uses $z=[0.25,0.50,0.85,1.20,1.60,2.00]$ m and the sparser three-sensor layout uses $z=[0.5,1.25,2.0]$ m. All three layouts include an outlet-temperature measurement at $z=2.0$ m, while outlet conversion is included as a separate scalar measurement. Measurement noise is added as

$$\begin{array}{cc}\hfill T_{meas}(z_i,t_k)& =T(z_i,t_k)+{ϵ}_{T,i,k},{ϵ}_{T,i,k}\sim \mathcal{N}(0,{\sigma }_T^2),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill X_{meas}\left(t_k\right)& =X_{out}\left(t_k\right)+{ϵ}_{X,k},{ϵ}_{X,k}\sim \mathcal{N}(0,{\sigma }_X^2),\hfill \end{array}$$

(16)

where ${\sigma }_T$ and ${\sigma }_X$ are the standard deviations of the temperature and conversion measurement noise, respectively; their numerical values for the base and stress-test cases are defined in Section 2.5.

2.3.2. Objective Function

At each time step, the unknown activity profile is represented by a vector of activity nodes,

$${\widehat{\mathit{a}}}_k={[{\widehat{a}}_{1,k},{\widehat{a}}_{2,k},\dots ,{\widehat{a}}_{m,k}]}^T.$$

(17)

The activity profile is estimated by solving a constrained nonlinear least-squares problem:

$$\begin{array}{cc}\hfill {\Phi }_k\left(\widehat{\mathit{a}}\right)& ={∥{\displaystyle \frac{{\mathit{T}}_{pred}\left(\widehat{\mathit{a}}\right)-{\mathit{T}}_{meas,k}}{{\sigma }_T}}∥}_2^2+{∥{\displaystyle \frac{X_{pred}\left(\widehat{\mathit{a}}\right)-X_{meas,k}}{{\sigma }_X}}∥}_2^2\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{\lambda }_s{∥\mathit{D}\widehat{\mathit{a}}∥}_2^2+{\lambda }_p{∥\widehat{\mathit{a}}-{\widehat{\mathit{a}}}_{k-1}∥}_2^2,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\widehat{\mathit{a}}}_k& =arg\underset{0.05\le \widehat{\mathit{a}}\le 1}{min}{\Phi }_k\left(\widehat{\mathit{a}}\right).\hfill \end{array}$$

(19)

Here $\mathit{D}\in {\mathbb{R}}^{(m-1)\times m}$ is the first-order finite-difference matrix that penalizes axial gradients in the estimated activity profile,

$$\mathit{D}=\left[\begin{array}{ccccc}-1& 1& 0& \dots & 0\\ 0& -1& 1& \dots & 0\\ ⋮& & \ddots & \ddots & ⋮\\ 0& \dots & 0& -1& 1\end{array}\right].$$

(20)

The first two terms fit the measured outputs; the third term discourages spatially oscillatory activity estimates; and the fourth anchors the trial profile to the previous estimate, reflecting the slow timescale of catalyst deactivation relative to the measurement sampling interval. The lower activity bound of 0.05 maintains numerical conditioning while allowing severe deactivation to be represented. The case-study weights are ${\lambda }_s=0.7$ and ${\lambda }_p=15$. These are algorithmic regularization parameters rather than covariance-derived quantities: ${\lambda }_s$ controls axial smoothness and ${\lambda }_p$ controls step-to-step persistence. Their contribution is examined through the ablation study in Section 4.5. The constrained least-squares problem is solved using SciPy’s Trust-Region Reflective algorithm with ${ϵ}_x={ϵ}_f={ϵ}_g={10}^{-8}$, a maximum of 250 function evaluations, and a warm start from the preceding estimate.

2.3.3. Estimator Performance Metrics

Estimator performance is evaluated using three scalar metrics. The activity-profile root mean square error,

$${\mathrm{RMSE}}_a=\sqrt{{\displaystyle \frac{1}{mN}}\sum _{k=0}^{N-1}\sum _{i=1}^m{\left({\widehat{a}}_{i,k}-a_{i,k}\right)}^2},$$

(21)

measures the overall accuracy of the reconstructed activity field across all nodes and time steps. The outlet-conversion prediction error,

$${\mathrm{RMSE}}_X=\sqrt{{\displaystyle \frac{1}{N}}\sum _{k=0}^{N-1}{\left({\widehat{X}}_{out,k}-X_{out,k}\right)}^2},$$

(22)

and the outlet-temperature prediction error,

$${\mathrm{RMSE}}_T=\sqrt{{\displaystyle \frac{1}{N}}\sum _{k=0}^{N-1}{\left({\widehat{T}}_{out,k}-T_{out,k}\right)}^2},$$

(23)

quantify how well the estimated activity profile reproduces the corresponding noise-free synthetic plant outputs. Conversion errors are reported in percentage points and temperature errors in kelvin.

2.3.4. Uniform-Activity Comparator

To test whether spatial temperature measurements add useful information beyond a bulk activity estimate, a conversion-only comparator is evaluated in every noise-and-sensor scenario. The comparator assumes uniform catalyst activity, ${\widehat{\mathit{a}}}_{unif,k}={\widehat{a}}_{unif,k}\mathbf{1}$, and estimates the scalar activity from outlet conversion:

$${\widehat{a}}_{unif,k}=arg\underset{0.05\le \widehat{a}\le 1}{min}{\left({\displaystyle \frac{X_{pred}\left(\widehat{a}\mathbf{1}\right)-X_{meas,k}}{{\sigma }_X}}\right)}^2+{\lambda }_p{\left(\widehat{a}-{\widehat{a}}_{unif,k-1}\right)}^2.$$

(24)

This comparator uses no axial temperature measurements and cannot reconstruct a spatial activity gradient. It therefore isolates the additional information provided by the proposed spatial estimator rather than serving as a no-estimation control.

2.4. Self-Optimizing Operating Layer

The self-optimizing layer updates reactor operating targets based on the estimated catalyst condition. For the heat-exchanged fixed-bed reactor in the present case study, the decision vector consists of the inlet gas temperature and the coolant inlet temperature:

$$\mathit{u}={[T_{in},T_{c,in}]}^T.$$

(25)

The optimization objective is to maximize phthalic anhydride yield while penalizing large target moves away from the nominal operating point:

$$J=Y_{PA}+w_S{S}_{PA,out}-{\lambda }_u{\parallel \mathit{u}-{\mathit{u}}_{ref}\parallel }_2^2,$$

(26)

where $Y_{PA}=X_{out}{S}_{PA,out}$ is the phthalic anhydride yield, $w_S>0$ is a small weight on outlet selectivity, and ${\lambda }_u>0$ penalizes movement away from the nominal reference targets. Although selectivity already enters implicitly through yield, the secondary term $w_S{S}_{PA,out}$ provides an additional optimization gradient at high-conversion conditions where yield is relatively flat but selectivity still varies; $w_S$ is kept small so that yield remains the dominant objective. The nominal fixed-target reference is ${\mathit{u}}_{ref}={[625,625]}^T$ K, which is also the baseline for fixed-target operation (Section 4.3). The hot-spot safety limit

$$T_{max}^{safe}=730\mathrm{K}$$

(27)

is enforced by excluding any grid-search candidate with $T_{max}>T_{max}^{safe}$ via a large infeasibility penalty, so the constraint is effectively treated as a hard bound within the search space. The case-study values are $w_S=0.01$ and ${\lambda }_u=2\times {10}^{-4}$ K⁻², and the thermal-target search covers $T_{in},T_{c,in}\in [620,632]$ K at 1-K resolution (169 candidate points total).

Self-optimized operation calculates updated values of $T_{in}$ and $T_{c,in}$ from the estimated activity profile using a two-dimensional grid search over the allowable thermal-target range:

$$[T_{in}^{*}\left(t\right),T_{c,in}^{*}\left(t\right)]=arg\underset{\mathit{u}\in \mathcal{U}}{max}J\left(\widehat{a}(z,t),T_{in},T_{c,in}\right),$$

(28)

where $\mathcal{U}$ denotes the feasible set of inlet temperature targets.

2.5. Performance and Robustness Evaluation Procedure

The primary performance endpoint is the optimization-value gap $\Delta Y_{PA}$, measuring the yield shortfall of the estimated-activity policy relative to the true-activity upper bound; secondary metrics are ${\mathrm{RMSE}}_a$, ${\mathrm{RMSE}}_X$, and ${\mathrm{RMSE}}_T$. The main estimation and target-update demonstration uses $N=15$ sampling steps. Each robustness scenario uses $N=10$ sampling steps; this shorter sequence reduces computational cost while retaining the full prescribed progression from fresh to final deactivated activity. Because no plant historian data are used, the noise levels are sensitivity-analysis assumptions rather than plant-specific instrument specifications. The baseline case uses ${\sigma }_T=4$ K and ${\sigma }_X=1$ percentage point; lower- and higher-noise cases use $({\sigma }_T,{\sigma }_X)=(2\mathrm{K},0.5\mathrm{pp})$ and $(6\mathrm{K},1.5\mathrm{pp})$, respectively. The number of axial temperature sensors is varied across $s=3,4,6$. Each scenario is repeated using 10 independent Gaussian-noise realizations (seeds 11, 22, …, 110), and the reported values are means and sample standard deviations across those trials.

Additional tests remove each regularization term individually (ablation study), compare performance against a conversion-only uniform-activity estimator, and perturb all plant kinetic pre-exponential factors by $-5\%$, $+5\%$, and $+10\%$ while retaining nominal kinetics in the estimator (model–plant mismatch study).

The optimization-value gap quantifies the performance loss attributable to estimation error:

$$\Delta Y_{PA}=Y_{PA}^{\mathrm{true}-\mathrm{opt}}-Y_{PA}^{\mathrm{est}-\mathrm{opt}},$$

(29)

where $Y_{PA}^{\mathrm{true}-\mathrm{opt}}$ is the phthalic anhydride yield obtained by optimizing with the true (known) activity profile, serving as the ideal upper bound, and $Y_{PA}^{\mathrm{est}-\mathrm{opt}}$ is the yield achieved by the practical digital-twin policy that uses only the estimated activity profile.

3. Case Study and Reactor Model

The proposed framework is demonstrated using the oxidation of o-xylene to phthalic anhydride (PA) in a heat-exchanged fixed-bed reactor. Phthalic anhydride is one of the most important aromatic intermediates in the chemical industry, with global production exceeding one million tonnes per year, primarily used in the manufacture of plasticizers, dyes, and resins [6,18]. The reaction system is strongly exothermic and has been widely used as a benchmark for fixed-bed catalytic reactor modeling, temperature-profile prediction, and safe-operation studies [5,6,19]. The desired reaction converts o-xylene to phthalic anhydride, while undesired reactions lead to over-oxidation and carbon-oxide formation. Reactor performance is strongly affected by temperature, catalyst activity, and heat removal.

3.1. Process Description

A dilute o-xylene feed is mixed with air, where oxygen is present in large excess and nitrogen is treated as an inert carrier. The gas mixture enters a bank of catalyst-filled tubes surrounded by molten-salt coolant. The catalyst is a vanadia–titania formulation commonly used for selective oxidation of o-xylene to phthalic anhydride. Heat generated by the strongly exothermic reaction network is removed through the tube wall to the coolant, and the reactor is operated under co-current coolant-flow conditions in the present model.

The process measurements assumed available to the digital twin are sparse axial gas temperatures and outlet conversion. These measurements are representative of the type of information used in model-based fixed-bed reactor studies and activity-estimation work [6,21]. The hidden state of interest is the axial catalyst activity profile $a(z,t)$, which modifies the apparent kinetic rates and is not measured directly. The self-optimizing layer uses the estimated activity profile to update the inlet gas temperature and coolant inlet temperature while respecting a hot-spot temperature constraint.

Figure 2 shows the reduced reaction network used in the digital-twin model. The network is intentionally simpler than the full industrial oxidation chemistry, but it preserves the key reaction-engineering tradeoff emphasized in previous o-xylene oxidation studies: the desired partial oxidation pathway competes with over-oxidation routes that reduce selectivity and increase heat release [6,19].

The simplified reaction network used in this study consists of three reactions:

$$\begin{array}{cc}\hfill A+3{\mathrm{O}}_2& \to B+3{\mathrm{H}}_2\mathrm{O},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill B+7.5{\mathrm{O}}_2& \to 8{\mathrm{CO}}_2+2{\mathrm{H}}_2\mathrm{O},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill A+10.5{\mathrm{O}}_2& \to 8{\mathrm{CO}}_2+5{\mathrm{H}}_2\mathrm{O},\hfill \end{array}$$

(32)

where A denotes o-xylene and B denotes phthalic anhydride. Reaction (30) is the desired partial oxidation; reactions (31) and (32) represent undesired over-oxidation of PA and direct combustion of o-xylene, respectively.

3.2. Reactor Configuration

The case-study reactor is an industrial multi-tube fixed-bed reactor for the partial oxidation of o-xylene to phthalic anhydride over a $V_2{O}_5/Ti{O}_2$ catalyst [4,23]. The reactor tubes are surrounded by molten-salt coolant, which removes the heat generated by the strongly exothermic oxidation reactions. The model represents one representative tube, while coolant heat transfer is scaled using the total number of tubes. The parameter values used in the simulation are summarized in Table 2; they follow the benchmark PA reactor model and kinetic data reported in the literature [6,19].

The reactor is modeled as a one-dimensional pseudo-homogeneous plug-flow reactor. The axial coordinate is denoted by z, with $z=0$ at the reactor inlet and $z=L$ at the reactor outlet. The ordinary differential equation (ODE) state vector integrated along the reactor is

$$\mathit{s}\left(z\right)={[F_A,F_B,F_C,F_N,T,T_c,P]}^T,$$

(33)

where $F_A$, $F_B$, $F_C$, and $F_N$ are the molar flow rates of o-xylene, phthalic anhydride, carbon oxides, and inert nitrogen, respectively; T is the gas-phase temperature; $T_c$ is the coolant temperature; and P is the total pressure. The symbol $\mathit{s}$ is used here to distinguish the case-study ODE reactor state from the general framework state vector $\mathit{x}$ introduced in Section 2.1.3.

The total inlet molar flow rate per representative tube is calculated from the superficial velocity and tube cross-sectional area:

$$F_{T,0}={\displaystyle \frac{u_s{\rho }_g{A}_{tube}}{M_g}},$$

(34)

where $u_s$ (m h⁻¹) is the superficial gas velocity at inlet conditions and $A_{tube}=\pi D_t^2/4$ is the tube cross-sectional area. The inlet o-xylene and inert molar flow rates are initialized as

$$\begin{array}{cc}\hfill F_{A,0}& ={\displaystyle \frac{P_{A,0}}{P_0}}{F}_{T,0},\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill F_{N,0}& =\left(1-{\displaystyle \frac{P_{A,0}}{P_0}}\right)F_{T,0},\hfill \end{array}$$

(36)

with initial conditions $F_B\left(0\right)=F_C\left(0\right)=0$, $T\left(0\right)=T_{in}$, $T_c\left(0\right)=T_{c,in}$, and $P\left(0\right)=P_0$.

3.3. Kinetic Model

The rate constants follow the de Lasa kinetic form [6,19]:

$$\begin{array}{cc}\hfill k_{21}\left(T\right)& =21.07exp\left(10.59-{\displaystyle \frac{13587.68}{T}}\right),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill k_{22}\left(T\right)& =21.07exp\left(11.62-{\displaystyle \frac{15801.97}{T}}\right),\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill k_{23}\left(T\right)& =21.07exp\left(9.73-{\displaystyle \frac{14392.88}{T}}\right),\hfill \end{array}$$

(39)

where the subscript notation $k_{2i}$ follows de Lasa [19]: the first digit (2) identifies the reaction family (selective oxidation of o-xylene), and the second digit ($i=1,2,3$) indexes the individual reactions corresponding to (30)–(32), respectively. All rate constants have units of kmol/(kg_cat h kPa). Catalyst activity modifies the apparent rate for each reaction directly through the activity factor $a(z,t)$, so no separate effective rate-constant definition is required; the activity enters the rate expressions shown below.

The local organic partial pressures are computed from molar flows:

$$\begin{array}{cc}\hfill F_T& =F_A+F_B+F_C+F_N,\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill P_j& ={\displaystyle \frac{F_j}{F_T}}P.\hfill \end{array}$$

(41)

The activity-modified reaction-rate expressions per unit reactor length are

$$\begin{array}{cc}\hfill r_A& =-{\rho }_{L}a(z,t)\left[k_{21}\left(T\right)+k_{23}\left(T\right)\right]P_A,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill r_B& ={\rho }_{L}a(z,t)\left[k_{21}\left(T\right)P_A-k_{22}\left(T\right)P_B\right],\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill r_C& =8{\rho }_{L}a(z,t)\left[k_{23}\left(T\right)P_A+k_{22}\left(T\right)P_B\right],\hfill \end{array}$$

(44)

where ${\rho }_L={\rho }_b{A}_{tube}$ is the linear catalyst density.

3.4. Mole, Energy, Coolant, and Pressure Balances

The reactor mole balances are

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{F}_A}{dz}}& =r_A,\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{F}_B}{dz}}& =r_B,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{F}_C}{dz}}& =r_C,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{F}_N}{dz}}& =0.\hfill \end{array}$$

(48)

The o-xylene conversion is calculated as

$$X_A\left(z\right)={\displaystyle \frac{F_{A,0}-F_A\left(z\right)}{F_{A,0}}},$$

(49)

and the phthalic anhydride selectivity is calculated as

$$S_{PA}\left(z\right)={\displaystyle \frac{F_B\left(z\right)}{F_{A,0}-F_A\left(z\right)}}.$$

(50)

The gas-phase energy balance includes heat generation by the exothermic reactions and heat removal to the molten salt coolant:

$${\displaystyle \frac{dT}{dz}}={\displaystyle \frac{Q_{gen}-Q_{cool}}{u{\rho }_g{C}_{p,g}}},$$

(51)

where

$$Q_{gen}={\rho }_{b}a(z,t)\left[k_{21}{P}_A|\Delta H_{21}|+k_{22}{P}_B|\Delta H_{22}|+k_{23}{P}_A|\Delta H_{23}|\right]$$

(52)

and

$$Q_{cool}={\displaystyle \frac{4U(T-T_c)}{D_t}}.$$

(53)

The heats of reaction are $\Delta H_{21}=-307,122$ kcal/kmol, $\Delta H_{22}=-783,700$ kcal/kmol, and $\Delta H_{23}=-1,090,822$ kcal/kmol [19]. The value of $\Delta H_{22}$ is confirmed from Hess’s law: $\Delta H_{22}=\Delta H_{23}-\Delta H_{21}=(-1,090,822)-(-307,122)=-783,700$ kcal/kmol, consistent with the stoichiometric relationship between reactions (30)–(32).

For co-current operation, the coolant energy balance is

$${\displaystyle \frac{d{T}_c}{dz}}={\displaystyle \frac{t_n\pi D_{t}U(T-T_c)}{w_c{C}_{p,c}}}.$$

(54)

The pressure drop is represented using an Ergun-type expression [24]:

$${\displaystyle \frac{dP}{dz}}=-{\displaystyle \frac{G^2(1-ϵ)}{{\rho }_g{D}_p{ϵ}^3}}\left[{\displaystyle \frac{150(1-ϵ)}{R{e}_p}}+1.75\right].$$

(55)

3.5. Model Assumptions and Validation

The main assumptions are: one-dimensional plug flow, negligible axial dispersion and radial gradients, pseudo-homogeneous catalyst bed, oxygen in large excess (so that oxygen partial pressure is approximately constant and does not appear explicitly in the rate expressions), constant gas properties, slow catalyst deactivation relative to the reactor residence time, and quasi-steady reactor behavior for any fixed activity profile. The timescale separation between deactivation and residence time is well satisfied in the present case: the reactor residence time is on the order of seconds, whereas the deactivation sequence used for estimation spans tens of time steps representing hours to days of operation, so the quasi-steady assumption is appropriate.

3.5.1. Forward-Model Validation

The model was validated under benchmark conditions with $T_{in}=T_{c,in}=625$ K, $P_{A,0}=0.9322$ kPa, and co-current coolant flow, matching the conditions reported in [6,19]. The model predictions for maximum temperature, outlet conversion, and PA selectivity were compared against the literature benchmark values; agreement was satisfactory across all three outputs. Table 3 summarizes the validation results.

3.6. Catalyst Activity-Profile Representation

Catalyst deactivation in fixed-bed reactors is often spatially nonuniform, with activity loss concentrated near the inlet or outlet depending on the deactivation mechanism. To capture spatial variation while keeping the estimation problem tractable, the continuous activity field is approximated by piecewise-linear interpolation between a finite number of axial activity nodes:

$$a(z,t)\approx {\mathcal{I}}_{\mathrm{lin}}\left[a_1\left(t\right),a_2\left(t\right),\dots ,a_m\left(t\right)\right],$$

(56)

where m is the number of activity nodes, ${\mathcal{I}}_{\mathrm{lin}}$ denotes piecewise-linear interpolation along the reactor axis, and the nodes are placed at uniformly spaced axial positions $z_i=(i-1)L/(m-1)$ for $i=1,\dots ,m$. In the present case study $m=4$, giving nodes at $z=[0,0.67,1.33,2.0]$ m; this choice provides sufficient spatial resolution to represent inlet-to-outlet activity gradients while keeping the estimation problem tractable given the five available measurements ($s=4$ temperatures plus one conversion). The activity vector is

$$\mathit{a}\left(t\right)={[a_1\left(t\right),a_2\left(t\right),\dots ,a_m\left(t\right)]}^T,$$

(57)

with estimation bounds

$$0.05\le a_i\left(t\right)\le 1.0,$$

(58)

where the lower bound of 0.05 prevents numerical singularity while permitting severe local deactivation to be represented, consistent with Section 2.3.

To characterize the effect of activity loss on reactor behavior and to verify observability before running the estimator, four representative axial activity-profile cases are examined: (i) fresh catalyst ($a_i=1$ everywhere), (ii) uniformly deactivated catalyst, (iii) linear outlet-side deactivation, and (iv) localized outlet-side activity loss. The forward-model responses for these cases are shown in Figure 3 and Figure 4. For the main estimation and self-optimization study, a synthetic time-evolving deactivation sequence of $N=15$ steps is used. The activity profile advances monotonically from fresh catalyst ($\mathit{a}={[1,1,1,1]}^T$) to a final deactivated state ${\mathit{a}}_{final}={[0.90,0.75,0.60,0.45]}^T$ (stronger deactivation toward the outlet, representing cumulative thermal and product exposure), with a nonlinear severity progression: at step k, ${\mathit{a}}_k=\mathbf{1}+{(k/14)}^{1.35}({\mathit{a}}_{final}-\mathbf{1})$. This synthetic trajectory provides a controlled severe-deactivation scenario for evaluating the estimator across a broad activity range; it does not represent a calibrated physical deactivation law.

4. Results and Discussion

All results in this section are obtained from a sequential simulation study using synthetic measurements generated from the benchmark forward model. The nominal estimator and target-update calculations use the same model structure, while additional kinetic-mismatch tests perturb the synthetic plant to probe sensitivity to model–plant error. The reported estimation and optimization figures should therefore be interpreted as simulation-based proof-of-concept results rather than plant-validated performance; factors that would reduce performance in deployment are discussed in Section 4.9.

This section presents and discusses the simulation results in six parts. Section 4.1 reports the forward-model validation and confirms that catalyst deactivation is observable from the available measurements. Section 4.2 presents the activity-profile estimation results. Section 4.3 demonstrates the self-optimizing operation and compares fixed-target, true-activity, and estimated-activity policies. Section 4.4 examines robustness to measurement noise and sensor density. Section 4.8 discusses novelty positioning, and Section 4.9 discusses limitations and future work.

4.1. Forward-Model Validation and Activity Observability

Under benchmark fresh-catalyst conditions, the model predicted $T_{max}=676.9$ K, $X_{out}=88.4\%$, and $S_{PA}=77.9\%$ (Table 3). These values are consistent with the benchmark results reported in the literature [6,19], confirming that the model is suitable as the physics engine of the digital twin. The temperature profile increases monotonically from inlet to exit under the present configuration (dilute feed, co-current molten-salt coolant, $T_{in}=T_{c,in}=625$ K), so the maximum temperature coincides with the reactor outlet rather than a mid-bed hot spot; this behavior is consistent with the benchmark model of de Lasa [19] at these operating conditions. Figure 5 shows the predicted axial temperature profile, confirming the monotonic rise from inlet to outlet and the absence of a mid-bed hot spot under these benchmark conditions.

Uniform catalyst deactivation was simulated at $a\in \{1.0,0.8,0.6,0.4\}$. As activity decreased, both the reactor temperature rise and outlet conversion decreased systematically, as shown in Figure 3. At $a=0.4$, the model predicted $T_{max}=639.8$ K and $X_{out}=26.9\%$, corresponding to a maximum-temperature reduction of approximately 37 K and a conversion reduction of approximately 61.5 percentage points relative to the fresh-catalyst case (Table 4). These results confirm that catalyst deactivation produces large, monotonic, and readily distinguishable changes in both thermal and compositional outputs across all four deactivation levels, establishing that the hidden activity state is strongly observable from the available process measurements. Figure 4 further shows that the representative axial activity profiles considered—uniform, linear outlet-side, and localized outlet-side deactivation—produce distinct axial temperature signatures, confirming that spatial activity gradients along the bed are also recoverable from the sparse temperature measurements used in the estimator.

4.2. Activity-Profile Estimation Results

The activity estimator was evaluated across 10 independent noise realisations (seeds 11–110) using a 15-step deactivation sequence under the baseline noise level (${\sigma }_T=4$ K, ${\sigma }_X=1$ percentage point). The estimator represented the activity profile using four axial nodes at $z=[0,0.67,1.33,2.0]$ m. The final true activity profile was

$${\mathit{a}}_{true}=[0.90,0.75,0.60,0.45],$$

(59)

representing stronger outlet-side deactivation consistent with greater thermal exposure near the bed exit. The mean final estimated profile across the 10 seeds was $\widehat{\mathit{a}}\approx [0.874,0.773,0.543,0.531]$. The full estimation performance metrics are reported in Table 5. The baseline activity-profile RMSE averaged across all 10 seeds and 15 steps was 0.075; the ${\mathrm{RMSE}}_X$ was 0.99 percentage points and the ${\mathrm{RMSE}}_T$ was 1.85 K, confirming that the estimated activity profile reproduces both outlet measurements with good fidelity despite the sparse sensor configuration. These results are important because the optimizer never receives the true catalyst activity profile; it receives only the estimated state reconstructed from four axial temperature measurements and one outlet-conversion measurement. In this sense, the estimator converts routinely available process measurements into actionable catalyst-state information, avoiding reliance on direct catalyst characterization during normal operation. The ensemble-mean profile snapshots in Figure 6 show the progressive deepening of the axial deactivation gradient: by step 14, the mean estimated profile declines from $\widehat{a}\approx 0.87$ at the inlet to $\widehat{a}\approx 0.53$ at the outlet. Thus, the estimator captures the main inlet-to-outlet deactivation trend, but the pointwise activity estimates remain noisy and partially biased because only sparse axial temperatures and one outlet-conversion measurement constrain the four-node activity profile. The measured and predicted outlet temperature and conversion fits are shown in Figure 7; both outputs track the simulated measurements closely with no systematic bias. The final smoothed profile in Figure 8 provides the operating-relevant activity distribution passed to the optimizer.

4.3. Self-Optimizing Operation Results

Under fixed operation, the reactor was maintained at $T_{in}=T_{c,in}=625$ K throughout the deactivation sequence. As catalyst activity declined to the final profile $[0.90,0.75,0.60,0.45]$, fixed-target operation progressively lost performance; the final phthalic anhydride yield under fixed targets dropped to approximately 46.3%.

The self-optimizing control results are from an illustrative single run (seed 44). Using the estimated activity profile, the optimizer updated the thermal targets to $T_{in}^{*}=632$ K and $T_{c,in}^{*}=632$ K over the deactivation sequence, recovering the final PA yield to approximately 62.1% while maintaining the maximum reactor temperature below the 730 K safety limit. The overall performance outcomes are summarized in Table 6.

Table 7 compares three policies at the final deactivated state under the baseline noisy measurement scenario (means over 10 noise realisations): (i) fixed-target operation ($\approx 46.3\%$), (ii) ideal optimization using the true activity profile (upper bound, $\approx 62.1\%$), and (iii) practical optimization using the estimated activity profile ($\approx 61.9\%$). The estimated-activity optimization falls within $0.14\pm 0.20$ percentage points of the true-activity optimum on average. Expressed as recovery efficiency, the true-activity optimum provides 15.75 percentage points of recoverable yield improvement over fixed-target operation, while the estimated-activity policy recovers 15.61 percentage points, or approximately 99.1% of the available performance recovery. This indicates that estimation uncertainty has little practical effect on the optimized operating target in the baseline scenario.

Figure 9. Self-optimized inlet gas and coolant temperature targets calculated from the estimated activity profile.

Figure 9. Self-optimized inlet gas and coolant temperature targets calculated from the estimated activity profile.

Figure 10. Performance recovery under estimated-activity self-optimization compared with fixed-target operation.

Figure 10. Performance recovery under estimated-activity self-optimization compared with fixed-target operation.

Figure 11. Maximum reactor temperature under fixed-target and self-optimized operation. The shaded blue band shows the temperature margin between the two policies; the red shaded zone above the dashed line marks the unsafe region ($T_{max}>730$ K). The self-optimized policy remains below the safety limit throughout the deactivation sequence.

Figure 11. Maximum reactor temperature under fixed-target and self-optimized operation. The shaded blue band shows the temperature margin between the two policies; the red shaded zone above the dashed line marks the unsafe region ($T_{max}>730$ K). The self-optimized policy remains below the safety limit throughout the deactivation sequence.

To connect the self-optimizing target update with safe-operation analysis, Figure 12 shows the two-dimensional operating map of the final deactivated reactor state as a function of gas and coolant inlet temperatures. The color field represents the maximum reactor temperature predicted by the model; white contours mark lines of constant PA yield (computed from the model at 2 percentage-point intervals); and the dashed black contour marks the 730 K hot-spot safety boundary. The step-by-step self-optimizing control trajectory is overlaid as a blue line, starting from the fixed nominal point (625 K, 625 K) and ending at the final self-optimized target. The map confirms that the trajectory moves into a region of substantially higher PA yield while remaining safely below the 730 K constraint throughout. This visualization follows the operating-window logic used in parametric sensitivity studies of o-xylene oxidation [6], but here it is generated from the deactivated catalyst state as estimated by the digital twin rather than from the fresh-catalyst model.

4.4. Robustness to Measurement Noise and Sensor Density

The stress tests showed that estimator and optimizer performance depend systematically on both measurement quality and sensor density (Table 8; Figure 13). All scenarios were evaluated over 10 sampling steps (seeds 11–110) using the estimator settings (${\lambda }_s=0.7$, ${\lambda }_p=15$) and optimizer settings ($w_S=0.01$, ${\lambda }_u=2\times {10}^{-4}$ K⁻²) from the main demo; reported values are means ± standard deviations across 10 independent noise realisations. Regarding noise sensitivity: in the low-noise case (${\sigma }_T=2$ K, ${\sigma }_X=0.5$ pp), the ${\mathrm{RMSE}}_a$ was $0.055\pm 0.009$ and the yield gap was $0.057\pm 0.12$ percentage points; in the baseline case (${\sigma }_T=4$ K, ${\sigma }_X=1$ pp), the yield gap was $0.14\pm 0.20$ percentage points; and in the high-noise case (${\sigma }_T=6$ K, ${\sigma }_X=1.5$ pp), the yield gap was $0.11\pm 0.20$ percentage points. The self-optimizing benefit over fixed-target operation is preserved across all three noise scenarios.

Regarding sensor density: the dense six-sensor configuration ($s=6$) achieved ${\mathrm{RMSE}}_a=0.070\pm 0.011$ with a yield gap of $0.11\pm 0.15$ percentage points, closely approaching the ideal true-activity optimum. The sparse three-sensor case ($s=3$) yielded ${\mathrm{RMSE}}_a=0.080\pm 0.012$ and a yield gap of $0.17\pm 0.24$ percentage points. The estimated-activity policy remained superior to fixed-target operation across all tested noise levels and sensor configurations.

4.5. Regularisation Ablation Study

To isolate the contribution of each regularisation term in the estimator objective, three ablation scenarios were evaluated at the baseline sensor configuration (4 axial temperature sensors, ${\sigma }_T=4$ K): (i) no spatial smoothing (${\lambda }_s=0$, ${\lambda }_p=15$), (ii) no temporal prior (${\lambda }_s=0.7$, ${\lambda }_p=0$), and (iii) no regularisation (${\lambda }_s=0$, ${\lambda }_p=0$). All 10 noise seeds were used per scenario. Table 9 reports the results.

The ablation results demonstrate a strong contribution from the temporal prior (${\lambda }_p=15$): removing it raises the activity-profile RMSE from 0.075 to 0.114 (a 52% increase), while removing only spatial smoothing causes a minor increase to 0.077. Removing both terms gives RMSE $=0.120$ (60% above the full-regularisation baseline). The temporal prior is the dominant regulariser because the per-step signal-to-noise ratio is below 0.4; without it, the optimizer fits noise directly, producing large step-to-step jumps. Spatial smoothing provides a smaller benefit because the piecewise-linear profile interpolation already partially couples adjacent nodes. The yield gap is essentially unchanged across all four configurations, consistent with the finding in Section 4.6 that the 1 K grid optimizer does not translate incremental accuracy improvements into measurable yield gains.

4.6. Comparison with Uniform-Activity Baseline

Table 10 compares the moving-window estimator against the uniform-activity baseline in every noise-and-sensor scenario. The uniform baseline uses only the outlet-conversion measurement and assumes spatially uniform deactivation; the moving-window estimator additionally incorporates sparse axial temperature measurements to resolve the spatial activity profile.

The activity RMSE of the uniform baseline was approximately 0.092 across all noise-and-sensor scenarios, reflecting the fact that uniform deactivation cannot be distinguished from the mean deactivation level of the true profile using a single outlet measurement. The moving-window estimator achieved lower activity RMSE than the uniform baseline across all tested scenarios, because the additional axial temperature measurements allowed it to resolve the spatial gradient more accurately than a single-output uniform estimator. Despite this, the moving-window estimator provides spatial information that the uniform baseline cannot: it identifies which section of the bed is most active, which is required for hot-spot safety monitoring and localized inlet-condition adjustments. In the optimization trials, the two-dimensional 1 K-resolution grid search was not fine enough to translate the spatial activity advantage into a consistent yield improvement; a finer-resolution or gradient-based optimizer would be expected to reveal larger gains from the spatial estimate.

4.7. Model-Mismatch Sensitivity

The nominal results reported above use a matched-model synthetic plant: the estimator and target optimizer use the same model structure that generates the nominal noisy measurements. To assess sensitivity to model–plant mismatch, three additional scenarios were run in which all reaction pre-exponential factors in the synthetic plant were perturbed by $-5\%$, $+5\%$, and $+10\%$ while the estimator and optimizer continued to use nominal kinetics. Table 11 summarises the results.

Under nominal conditions (0% perturbation), the yield gap was $0.14\pm 0.20$ pp, consistent with the stress-test baseline. With $+5\%$ perturbation, the gap was $0.16\pm 0.22$ pp, essentially unchanged. With $+10\%$ perturbation, the gap increased to $0.45\pm 0.33$ pp; the higher absolute yield levels in this scenario ($\approx 67.7\%$ vs. $61.9\%$ nominal) reflect the more reactive perturbed plant rather than improved estimation. The $-5\%$ perturbation produced a gap of $0.11\pm 0.20$ pp, slightly below nominal. These results indicate that the estimated-activity policy remains within less than 0.5 pp of the true-activity optimum across the full tested mismatch range, evaluated on the same perturbed plant in each case. The result does not eliminate the need for plant validation, but it suggests that the target-update layer is relatively insensitive to moderate uniform kinetic pre-exponential-factor mismatch; its value comes from repeatedly adapting operating targets using the inferred catalyst state rather than from a single fixed nominal design calculation.

4.8. Novelty Positioning Against Existing Literature

As summarized in Table 1, catalyst activity estimation in fixed-bed reactors has been studied previously, including studies where temperature measurements are used to infer spatial activity profiles [21]. Reactor optimization and model-based digital design have also been reported for catalytic reactors, including the PA system [5,22]. The novelty of the present work is therefore not the introduction of a new kinetic model, nor the isolated concept of catalyst activity estimation. Instead, the contribution lies in demonstrating, in a simulation-based proof of concept, how an estimated catalyst activity profile is propagated directly into a self-optimizing target-update layer—closing the loop from sparse measurements to hidden-state estimation to updated operating targets. This distinction is important because measurements alone reveal that reactor performance has deteriorated, but they do not identify how the operating targets should be adjusted; the estimated catalyst state supplies the physical information needed to select catalyst-specific temperature targets.

In this framework, the estimated catalyst activity serves as the operative input to the optimization problem at each time step, not merely as a diagnostic output. This transforms the digital twin from a passive monitoring tool into an adaptive decision-support system that adjusts operating targets as the catalyst degrades. From an industrial perspective, maintaining operation near the catalyst-specific optimum could help delay conservative replacement decisions, reduce the economic penalty of progressive deactivation, and preserve yield without requiring extensive sensor retrofitting. The present work is scoped as a computational methodology contribution using a validated benchmark reactor model and synthetic measurements; plant deployment would require further validation against experimental or industrial data, as discussed in Section 4.9.

4.9. Limitations and Future Work

The five limitations listed below are ranked by their effect on confidence in the reported performance numbers. The most important is the simulation-only validation basis (Limitation 1), which means the reported RMSE and yield figures cannot be directly extrapolated to real plant conditions. The second most important is the simplified 1D model (Limitation 2), which constrains the practical scope. The remaining three (Limitation 3: coarse activity discretisation; Limitation 4: simplified estimator; Limitation 5: simplified optimizer) constrain the deployment pathway but do not affect the internal consistency of the simulation results.

Limitation 1: Simulation-only validation. The study is a simulation-based proof of concept. The nominal synthetic measurements are generated from the same benchmark model structure used by the estimator, with additive Gaussian noise superimposed, and the kinetic-mismatch study perturbs only the reaction pre-exponential factors. Therefore, the results do not cover all forms of structural model error that would be present in plant deployment. Heat-transfer degradation, non-Gaussian sensor noise, sensor bias, catalyst-batch variability, and structural kinetic uncertainty could each reduce performance relative to the reported values. Plant data, independent experimental measurements, or pilot-scale tests would be required to quantify these effects and assess deployment feasibility.

Limitation 2: One-dimensional pseudo-homogeneous reactor model. The reactor model is a one-dimensional pseudo-homogeneous model. This level of model detail is appropriate for repeated estimation and optimization because it is computationally efficient and reproduces the benchmark reactor behavior used in this study. However, it neglects radial temperature gradients, intraparticle diffusion limitations, detailed catalyst morphology, and possible nonuniform coolant-side effects. A more detailed two-dimensional heterogeneous model would be required for industrial reactor design or safety certification.

Limitation 3: Low-dimensional activity discretisation. The catalyst activity profile is represented by a small number of axial nodes. This low-dimensional representation improves numerical conditioning and avoids overfitting sparse measurements, but it cannot resolve fine-scale deactivation patterns. The smoother shape of the estimated profile relative to the true profile (Section 4.2) is an estimator artifact—a consequence of spatial regularization and limited sensor coverage—rather than a physical feature of the catalyst state. The estimated profile should therefore be interpreted as an effective activity distribution that supports operating decisions, not as a direct microscopic measurement of the catalyst state. These choices limit geometric fidelity and prevent fine-scale deactivation reconstruction.

Limitation 4: Single-step moving-window estimator. The estimator is a moving-window constrained activity estimator with spatial smoothing and temporal regularization. It is deliberately described as a moving-window estimator rather than a full nonlinear moving-horizon estimation (MHE) implementation [14,25]: the current formulation estimates activity independently at each time step using the previous estimate as a regularization anchor, rather than optimizing over a sliding window of past measurements while enforcing dynamic process constraints simultaneously. A rigorous MHE formulation would improve estimation consistency and enable principled tuning of the horizon length.

Limitation 5: Simplified steady-state optimizer. The optimization layer is a simplified steady-state operating-target update, not a plant-ready economic model predictive controller or real-time optimizer [25,26,27]. Only two manipulated variables—inlet gas temperature and coolant inlet temperature—are optimized via a coarse two-dimensional grid search. Practical deployment would additionally require dynamic actuator models, ramp-rate limits, feed-flow constraints, coolant-system hydraulics, economic cost terms, catalyst-life penalty functions, and formally verified plant safety margins. Taken together, these five limitations define the scope of the present work as a simulation-based proof-of-concept computational framework.

Future work can be organized into two tiers. Near-term methodological extensions include: (i) replacing the current step-by-step estimator with a rigorous sliding-window nonlinear moving-horizon estimation (MHE) formulation to improve estimation consistency and enable principled horizon-length tuning—directly addressing Limitation 4; (ii) replacing the grid-search optimizer with a model predictive control (MPC) layer that handles dynamic actuator constraints and multiple manipulated variables—addressing Limitation 5; and (iii) testing finer activity discretizations ($m>4$) to assess the resolution limit imposed by available sensor configurations. Longer-term priorities include: (iv) validation against pilot-scale or industrial plant data to assess robustness to structural model mismatch, sensor bias, and catalyst-batch variability—the primary validation gap identified in Limitation 1; and (v) extension of the optimization objective to include economic catalyst-life penalty terms and heat-transfer degradation, moving the framework toward full economic MPC for deactivating reactors.

5. Conclusions

Catalyst deactivation progressively degrades fixed-bed reactor performance, yet most reactor digital twins lack the ability to estimate the current catalyst condition and use it to adapt operating targets. This work addressed that gap by developing a catalyst-activity-aware self-optimizing digital twin framework for catalyst-deactivating fixed-bed reactors. The central idea is to treat catalyst activity as a hidden reactor state, estimate it sequentially from sparse process measurements, and use the estimated condition to update reactor operating targets rather than holding them fixed.

The framework was demonstrated using the oxidation of o-xylene to phthalic anhydride over a vanadia–titania catalyst in a heat-exchanged fixed-bed reactor. A one-dimensional pseudo-homogeneous reactor model was implemented and validated under benchmark fresh-catalyst conditions ($T_{in}=T_{c,in}=625$ K), reproducing the literature benchmark values of $T_{max}$, $X_{out}$, and $S_{PA}$ reported by de Lasa [19] and Zuluaga-Botero et al. [6]. Catalyst activity was incorporated as an axial activity profile modifying the apparent reaction rates, directly linking catalyst-state loss to outlet conversion, selectivity, heat release, and hot-spot location.

Simulation results confirmed that catalyst deactivation produces large, monotonic changes in axial temperature and outlet conversion. At 40% activity relative to fresh catalyst, $T_{max}$ decreased by approximately 37 K and $X_{out}$ by approximately 61.5 percentage points, establishing that the hidden activity state is strongly observable. The moving-window constrained estimator recovered the dominant inlet-to-outlet deactivation gradient with an activity-profile ${\mathrm{RMSE}}_a=0.075$ (mean over 10 seeds), an outlet-conversion ${\mathrm{RMSE}}_X=0.99$ percentage points, and an outlet-temperature ${\mathrm{RMSE}}_T=1.85$ K.

The estimated activity profile was supplied to a self-optimizing operating layer that updated the inlet gas and coolant temperature targets at each estimation step. Under the baseline noisy measurement scenario (${\sigma }_T=4$ K, ${\sigma }_X=1$ percentage point), the estimated-activity optimization raised the final phthalic anhydride yield from 46.3% under fixed-target operation to 61.9%, within 0.14 percentage points of the ideal upper bound obtained with the true activity profile, while maintaining the maximum reactor temperature below the 730 K safety limit. Relative to the true-activity optimum, the estimated-activity policy recovered approximately 99.1% of the yield improvement available above fixed-target operation, showing that the reconstructed catalyst state was sufficiently accurate for operating-target selection.

Robustness testing over noise levels ${\sigma }_T=2$–6 K and sensor configurations $s=3$–6 (each scenario evaluated over 10 noise realisations) showed that the estimated-activity policy remained superior to fixed-target operation across all tested scenarios, with yield gaps from the true-activity optimum ranging from $0.06$ percentage points (low-noise, four sensors) to $0.17$ percentage points (three sensors, baseline noise). A regularisation ablation study confirmed that the temporal prior (${\lambda }_p=15$) is the dominant regulariser, with its removal raising activity-profile RMSE by 52%. Kinetic mismatch testing showed that the estimated-activity policy remained within 0.5 pp of the true-activity optimum across the full tested range ($-5\%$ to $+10\%$ perturbation, evaluated on the perturbed plant), indicating robust performance across plausible levels of kinetic uncertainty.

All reported performance values were obtained from simulation. Although the kinetic-mismatch tests provide an initial sensitivity check, performance under real plant conditions, with structural model errors, sensor biases, non-Gaussian disturbances, and unmeasured heat-transfer changes, remains to be quantified through experimental or pilot-scale validation.

The main contribution of this work is the demonstration, in a simulation-based proof-of-concept framework, that constrained catalyst activity estimation and operating-target optimization can be integrated inside a unified digital twin loop to recover reactor performance under catalyst deactivation. The results suggest that catalyst-state-aware operation can convert sparse process measurements into practical operating decisions, with potential industrial value for preserving yield, delaying overly conservative catalyst replacement, and reducing the economic impact of progressive deactivation. Near-term future work should address full nonlinear moving-horizon estimation with a rigorous sliding-window formulation, model predictive control replacing the current grid-search optimizer, and extension to additional manipulated variables. Longer-term priorities include economic optimization with catalyst-life penalty terms, incorporation of heat-transfer degradation, and validation against pilot-scale or plant data to assess robustness to structural model mismatch.