A reliability-driven digital twin framework for maintenance scheduling: from asset criticality to hazard-based decision optimisation

2.1. Stage 1 – asset criticality ranking (Delphi-AHP-RPI)

The first stage determines which assets justify deeper modelling. The chosen method is a three-step procedure that integrates qualitative expert elicitation, a quantitative weighting model, and a percentile-based risk index. The fit is deliberate: criticality decisions in industrial maintenance are driven by both expert tacit knowledge and quantitative trade-offs, and no single method captures both. Fig. 2 summarises the three-step Delphi-AHP-RPI workflow.

Step 1a uses a Delphi study to elicit and rank candidate criteria. A panel of ten experts is recruited through purposive and snowball sampling, with each expert holding at least five years of relevant maintenance or asset management experience, in line with CREDES guidance recommending panel sizes between 10 and 18 [24]. Experts complete a semi-structured questionnaire that maps sub-criteria, for instance, Average Annual Failure Rate (AAFR), Downtime Per Failure (DPF), and Replacement Cost (RC), onto the nine main criteria proposed by Jasiulewicz-Kaczmarek et al. [25]. The decision rule is that Fleiss’s kappa is at least 0.50 across criteria, with up to three rounds of feedback before the most frequently selected criteria are accepted by a majority. A preliminary correlation matrix prunes redundant sub-criteria with correlations above an 80 % threshold to avoid double-counting.

Step 1b applies the Analytic Hierarchy Process (AHP) to convert expert judgements into priority vectors [26]. Pairwise comparisons use Saaty’s 1-9 intensity scale across a four-level hierarchy: goal, main criteria, sub-criteria, and asset alternatives. The number of pairwise judgements per matrix is $n\left(n-1\right)$/2, with n the number of compared elements. For each pairwise matrix, the priority vector is computed by column normalisation, row summation, and division by the number of columns; the maximum eigenvalue, ${\lambda }_{max}$, is obtained by multiplying the original matrix by the priority vector and dividing element-wise. Internal consistency is enforced through the Consistency Index and Consistency Ratio, with CR < 0.10 as the acceptance threshold. The geometric mean is used to aggregate the 10 experts’ ratings because it preserves the reciprocal nature of the Saaty scale:

Fig. 2Stage 1 asset criticality ranking. Delphi-AHP-RPI: expert elicitation

Step 1c integrates AHP weights with Failure Mode, Effects, and Criticality Analysis (FMECA) into a Risk Priority Index (RPI). Severity (S), Occurrence (O), and Detection (D) are scored on percentile-based 5-level rating scales drawn from the empirical distribution of each criterion in the asset population, so that thresholds reflect the data rather than analyst intuition. Occurrence is mapped to time-to-failure ranges; severity is mapped to criterion values, such as DPF (in hours) and RC (in monetary units). Detection is set to 1 when the framework is applied to assets already being monitored, which focuses the ranking on severity and occurrence. The Risk Priority Number and Risk Priority Index are:

where $P{V}_i$ is the priority vector at hierarchy level $i$. Tied assets are separated by the highest severity rating. The output of Stage 1 is a ranked asset list with explicit priority vectors at every hierarchy level, an audit trail of expert agreement (Fleiss’s kappa), and a CR < 0.10 check for every pairwise matrix. The asset at the top of the RPI ranking is carried into Stage 2.

2.2. Stage 2 – reliability structure modelling (BowTie hybrid)

Stage 2 explains how the failure logic of each critical asset is modelled using BowTie analysis. BowTie illustrates how threats lead to a top event (on the fault tree's left side) and how consequences develop (on the event tree’s right side). Fault Tree Analysis (FTA) determines the logical combinations of basic events that trigger the top event; BowTie adds preventive and mitigative barriers on top. Reliability Block Diagram (RBD) concepts organise barriers in series, parallel, or k-out-of-n arrangements, but BowTie applies a defence-in-depth principle: if any barrier succeeds, the process halts. This hybrid approach complies with IEC 61025 [22] and BS EN 61078 [23], offering a decision-focused perspective that connects causes, barriers, and outcomes. Fig. 3 summarises the three-step BowTie hybrid workflow.

Fig. 3Stage 2 reliability structure modelling workflow using a BowTie hybrid

Step 2a involves three steps. First, the top event is identified based on the critical failure mode identified in Stage 1, and a Fault Tree Analysis (FTA) is constructed. Threat failure rates ${\lambda }_i$ are derived from historical data.

Step 2b lists preventive barriers. Each barrier’s metrics, including probability of failure on demand (PFD_avg), failure per hour (PFH), human-error probability (HEP), or conditional probability ($p_j$), are transformed into a time-dependent reliability $R_j\left(t\right)$. For a barrier with a constant failure rate ${\lambda }_j$, this relationship is given by:

giving a failure probability:

For a constant error or conditional probability ($p_j$), the reliability is written as:

The probability that all preventive barriers on a threat pathway $i$ fail is written as:

Summing over threats yields the top-event frequency:

Step 2c models mitigative barriers in the same way. The consequence frequency is expressed as:

Two decision rules link Stage 2 to Stage 3. A single-point-of-failure (SPOF) rule flags any basic event whose unmitigated occurrence triggers the top event; such SPOFs are either reinforced with redundancy or escalated to Stage 3. A barrier-importance rule uses Risk Achievement Worth (RAW) and Barrier Importance Factor (BIF) to prioritise barriers, which are written as:

High-RAW or high-BIF barriers are reinforced or improved; low-impact barriers may be removed. The output is a list of dominant threats, critical components and high-impact barriers; these feed the hazard-based twin risk model in Stage 3.

2.3. Stage 3 – hazard-based twin risk model (Cox PH with residual learning)

Stage 3 estimates the probabilistic risk influencing the schedule. It uses a two-layer Cox proportional hazards (PH) model with residual-learning correction, structured in three stages: a fleet baseline, a condition-monitoring (CM) offset, and a cross-fitted machine-learning corrector. The model fitting is explicit. Cox PH is semi-parametric, manages right-censored lifetimes, and provides interpretable hazard ratios, which are important for transparency in maintenance decision-making [27], [28]. The two-layer split preserves industrial context (the full fleet) while learning from the smaller monitored subset where covariates are recorded. Residual learning addresses the standard limitation of linear Cox: nonlinear and interaction effects in covariates that the linear predictor cannot capture [28]. Fig. 4 summarises the three-step Cox PH and residual-learning workflow.

Step 3b expresses the base hazard for an asset $i$ with covariate vector $x_i$ is:

where $h_0\left(t\right)$ is the baseline hazard, and beta is the vector of log-hazard ratios. In the proposed framework, this is extended into a two-layer formulation. The fleet layer is fitted on $N$ units (e.g. 72 in the authors’ pilot dataset) with fleet-wide covariates $w_i$:

The CM layer is fitted on the monitored subset (typically a fraction of $N$) using $h_{0,F}\left(t\right)$ as a fixed offset, so that fleet-level information is preserved while the linear predictor learns from the monitored covariates $x_i$:

A residual corrector $z\left(x_i\right)$ is then trained out-of-fold (OOF) on the monitored subset in Step 3c. The OOF protocol partitions the monitored subset into $K$-folds, fits the CM-layer Cox on $K-1$ folds, generates residuals on the held-out fold, and trains a non-linear learner, for example, a gradient-boosted tree, to predict the residuals on the held-out fold. The final deployment model refits Cox jointly on beta and gamma, with $z\left(x_i\right)$ computed from the trained learner:

The baseline cumulative hazard $H_{0,F}\left(t\right)$ is estimated by the Breslow estimator:

where:

and $d_j$ is the number of events at time ( $t_j)$ and $Y_k\left(t_j\right)$ is the at-risk indicator. Covariates are $z$-standardised before fitting ($x_z=x-\mu /\sigma$); $\mu$ and $\sigma$ are stored as part of the deployment artefact, so that incoming sensor values are scaled consistently inside the DT.

Two assumption checks, in Step 3d, are mandatory before the model is exported. The proportional-hazards assumption is tested using Schoenfeld residuals [29] over time and by assessing the linearity of the log-cumulative hazard. Martingale residuals are inspected for systematic dependence on covariates [30]. If either check fails on a covariate, either a stratified-baseline Cox is fitted, or the residual-learning correction is enabled to absorb the non-linearity. The data leakage rule is firm: the residual target used to train z must always be out-of-fold, and the final $\gamma$ is refitted on the full monitored subset only after OOF predictions are in place.

Fig. 4Stage 3 hazard-based twin risk model workflow

Variable specification follows ISO 17359 indicators [21]. Condition-monitoring covariates include vibration RMS, temperature, torque, and oil debris or particle count; physics-informed derivations (engine torque from power and speed; corrected turbocharger speed using inlet-air temperature; inlet-air density from absolute pressure and temperature) are available as engineered covariates when sensor coverage permits. Maintenance covariates include maintenance count, time since last maintenance, event class (preventive, corrective, inspection), and observed MTTR. Operational covariates include load, rotational speed, ambient temperature, and operator shift. Censoring is treated as right-censored and assumed non-informative conditional on covariates; the lifecycle outcome is binary (failure vs censored), and the post-maintenance state is treated as as-good-as-new for the bearing-replacement scenario considered in the pilot.

The output of Stage 3 is a per-asset hazard $h_i\left(t\right)$, survival $S_i\left(t\right)$, risk multiplier ${\eta }_i=e^{\left(\beta x_i+\gamma z_i\right)}$, and a calibrated $H_{0,F}\left(t\right)$. These are the inputs to Stage 4.

2.4. Stage 4 – decision and scheduling optimisation (parametric simulation)

Stage 4 converts the hazard into decision-ready maintenance quantities, including rolling-window failure probability, expected downtime, cost, and availability. The first implementation is a parametric simulation over a two-dimensional decision grid (preventive-interval $\mathrm{\Delta }$ and risk-threshold $\theta$). For each candidate ($\mathrm{\Delta },\theta$), the simulation applies trigger logic, computes cost, downtime, failure-risk and availability outcomes, and identifies the best feasible policy under an availability constraint. This paper opts for a grid-based simulation as the first implementation because the decision space is low-dimensional and directly auditable. A Pareto front is then extracted to support a decision maker who must trade cost, downtime, and risk. Fig. 5 summarises the three-step decision and scheduling-optimisation workflow.

Fig. 5Stage 4 decision and scheduling optimisation workflow

In Step 4a, the decision quantities are derived from the hazard and survival functions of Stage 3. The probability of failure inside an interval $\mathrm{\Delta }$ for asset $i$ is:

Expected downtime, total cost, and availability are:

where $T$ is the analysis horizon (typically 365 days), $C_{PM}$ is the cost per preventive action, and $C_{CM}$ is the cost per corrective action. The optimisation is:

The decision space is searched over $\mathrm{\Delta }$ in [30, 365] days and $\theta$ between the 70th and 95th percentile of the fleet hazard distribution in Step 4b. The combined grid yields roughly 10³-10⁴ candidate policies. For each ($\mathrm{\Delta }$, $\theta$) pair in Step 4c, the simulation engine evaluates the asset over the analysis horizon: if $h_i\left(t\right)$ crosses $\theta$, an early preventive intervention is triggered; if the time since the last preventive action reaches $\mathrm{\Delta }$, a scheduled preventive action is triggered; otherwise, monitoring continues. Stochastic failure events are drawn from $S_i\left(t\right)$. The simulation aggregates total downtime, total cost, and availability per policy.

Pareto extraction, in Step 4d, can be applied after the grid simulation to identify the non-dominated candidate policies ($\mathrm{\Delta }$, $\theta$) across cost, availability, and expected number of failures. The recommended policy (${\mathrm{\Delta }}^{*},{\theta }^{*}$) is then selected by the asset manager based on risk appetite. Validation of the simulation engine is enforced through three checks: monotonicity (increasing $\mathrm{\Delta }$ within the tested grid drives $P_{fail}$ and downtime upward and availability downward), extreme values (very small $\mathrm{\Delta }$ must produce a PM cost spike; very large $\mathrm{\Delta }$ must produce a downtime spike), and analytical consistency checks against expected behaviour.

2.5. Stage 5 – DT integration and offline validation

The DT prototype is implemented as a four-layer stack, as outlined in Step 5a, with unidirectional data flow across the model, simulation, decision, and visualisation layers. The model layer hosts the hybrid Cox engine and the standardisation parameters. The simulation layer hosts the parametric optimisation. The decision layer applies the trigger logic and emits the recommendation. The visualisation layer presents hazard curves, survival curves, base versus corrected risk scores, and the recommended schedule with cost and downtime estimates. A dual-track display, based on Cox versus Cox with residual correction, is provided to support auditability: when the two paths disagree, the magnitude of $\gamma \cdot z\left(x\right)$ is reported to expose the contribution of the corrector. Fig. 6 summarises the three-step DT integration, calibration, and validation workflow.

The DT operates offline as an advisory tool. Implementing real-time streaming is planned for future development. Data flows only from the asset to the DT, with no feedback to plant control systems. This design choice maintains the framework’s auditability and supports operator-in-the-loop decision-making, while also isolating the analytical workflow from potential control-loop interactions.

Three validation pathways are used in Step 5b before deployment. Internal validation generates synthetic asset lifetimes from a known Cox model with controlled covariate perturbations ($\pm 1\sigma$, $\pm 2\sigma$); the predicted hazard and survival are scored against the synthetic ground truth using the Kullback-Leibler divergence [31]. Temporal validation trains the DT on years 1-4 of historical data and tests on year 5, scoring per-asset $P_{fail}\left({\mathrm{\Delta }}^{*}\right)$ using Harrell’s C-index, the Brier score, and the Expected Calibration Error (ECE) [32]-[34]. Sensitivity analysis perturbs the final $\beta$ by ±1 SE, perturbs final $\gamma$ similarly, and varies $C_{PM}$ and $C_{CM}$ by ±20 % to test the stability of the Pareto front and of the recommended (${\mathrm{\Delta }}^{*},{\theta }^{*}$) policy.

Fig. 6Stage 5 DT integration and offline validation workflow