Gearbox bearing crack growth prognostics and uncertainty quantification with physics-informed machine learning

This paper introduces the extreme theory of functional connections (X-TFC), a physics-informed machine learning algorithm, and tailors it to estimate the remaining useful life (RUL) of wind turbine gearbox bearings experiencing fatigue crack growth. Unlike purely data-driven methods, X-TFC embeds a physics model, based on Head's theory in this work, into its training objective. The core of X-TFC is a random-projection single-layer neural network trained via an extreme learning machine, which requires only limited damage progression data and solves for output weights with a least-squares optimization algorithm. A composite loss function balances the network's fit to observed degradation data against the residuals of the governing crack growth differential equation, ensuring the learned damage trajectory remains physically plausible. When applied to a vibration-based health-index (HI) dataset measured during the growth of a crack on the inner ring of a high-speed bearing in a wind turbine gearbox (Bechhoefer and Dubé, 2020), X-TFC achieves near-zero prediction bias. Even when trained on only the first 10 %–20 % of the damage progression data, with sufficient physics weighting its predictions remain monotonic and smooth, delivering high prognosability and trendability. To quantify the epistemic uncertainty, we employ a Monte Carlo ensemble of independently initialized X-TFC models trained on noise-perturbed data, which yields confidence intervals around each RUL estimate and captures both model-parameter and epistemic uncertainty. In addition to a vibration-based HI, we demonstrate that the proposed framework can be directly applied to a supervisory control and data acquisition (SCADA) data-based HI (Eftekhari Milani et al., 2026) measured during similar wind turbine gearbox bearing crack faults, preserving its accuracy and interpretability. This extension shows the versatility of our approach, which is applicable to bearings of multiple gearbox manufacturers, models, and ratings using only SCADA data. By integrating domain knowledge with machine learning, X-TFC offers a rapid, reliable tool for crack prognostics. Its adaptability to other bearing failure modes, such as pitch bearing ring cracks, positions X-TFC as a powerful enabler of data-driven, physics-informed asset management in the wind energy sector and beyond.

1 Introduction

The wind industry has witnessed tremendous growth over the past 2 decades, as reported by the Global Wind Energy Council, with the global cumulative installed capacity recently crossing the 1000 GW mark. Although the levelized cost of energy from wind and wind plant operation and maintenance (O&M) costs have fallen (Wiser et al., 2024), O&M still accounts for up to 35 % of the levelized cost of energy. Approximately half of wind plant O&M costs are from wind turbine O&M, a cost of USD 22 kW-year for land-based wind plants (Wiser et al., 2019). Turbine O&M costs, therefore, represent the single largest component of wind plant O&M costs and the primary source of potential O&M cost reductions. More reliable components and better O&M strategies, including for wind turbine blades and drivetrain components, have the potential to reduce premature failures, O&M costs, downtime, and thus the levelized cost of energy by up to 7 % by 2035 (Stehly et al., 2024). Additionally, an accurate reliability forecast provides crucial information for reducing O&M costs through design improvements, optimized operation strategies, and enhanced budgeting. Optimized operation strategies often include condition-based maintenance and prognostics, which rely on the prediction of the remaining useful life (RUL) of components to provide an early warning of failure and thereby allow maintenance to be scheduled proactively (Bechhoefer et al., 2021).

One challenging reliability problem that remains for wind turbine main bearings and bearings in gearboxes is axial cracking of the inner rings, examples of which are shown in Fig. 1a and b (Greco et al., 2022). In gearboxes, this failure mode in intermediate- and high-speed stage bearings accounts for the majority of uptower bearing replacements (Haus et al., 2024). “Axial” describes the orientation of these cracks as they align with the axis of rotation of the shaft on which the bearing is mounted. This type of crack is also commonly called a white-etching crack (WEC) due to the appearance of the steel microstructure when the cracked bearing cross-sections are polished, etched with chemicals, and examined under reflected light. These cracks tend to propagate to spalls or lead to a complete splitting of the bearing inner ring in gearboxes in as little as 5 % to 20 % of the bearing rating life (Greco et al., 2013). The conditions leading to WECs, the reasons for their apparent prevalence in wind turbine main bearings and gearboxes (Greco et al., 2022; Demas et al., 2023), and the effectiveness of countermeasures (Jensen et al., 2021) are now better understood through over a decade of research; however, accurate methods to predict the RUL of bearings with WECs still need to be developed. Early efforts focused on the development of a physics-based probability of failure model using frictional energy accumulation, which has been partially validated on a limited set of damaged bearings from a wind plant; however, the probability of failure and RUL were predicted only for the population of gearboxes rather than individual ones (Guo et al., 2020). The model was also used to explore the influence of bearing clearance and lubricant temperature (Clark et al., 2023). Another challenging reliability problem that has recently become more prevalent is fatigue cracking of pitch bearing outer rings in larger wind turbines, as shown in Fig. 1c, resulting in a replacement rate of several percent per year and a 10 % replacement rate for the population in just 5 years (Haus et al., 2025). These cracks can occur at bolt holes, fill holes, or locations of significant stiffness changes (Shapiro, 2017) and risk the blade detaching from the turbine. A typical pitch bearing replacement requires the use of a large, expensive crane and removal of the blade(s) or rotor, with a significant amount of downtime. However, research into the root causes of pitch bearing ring cracking has only just begun (Schleich, 2025). In rotating machinery, fatigue cracks in bearings, gears, or shafts are a common failure mode. Classical fracture mechanics models, such as linear elastic fracture mechanics, Head's theory, and dislocation theory, describe such crack growth under cyclic loading. However, no single crack growth law fits all phases of damage: Bechhoefer and Dubé (2020) showed that three high-cycle fatigue modes were needed to cover the primary crack propagation period. They found that a hybrid combination of modes was more accurate than any individual model.

Figure 1 Examples of bearing cracks and damage in wind turbines. (a) Gearbox bearing WEC and spalling (Greco et al., 2022). (b) Main bearing WEC and spalling (Greco et al., 2022). (c) Pitch bearing ring crack (Shapiro, 2017).

In recent years, new computational techniques have been developed to improve evaluation and prediction tasks in various fields of science. Among these, physics-informed machine learning (PIML) (Karniadakis et al., 2021) has emerged as a powerful paradigm for integrating data-driven modeling with known physical laws, often expressed using differential equations. By embedding governing equations and domain constraints directly into learning objectives, PIML methods improve robustness, data efficiency, and interpretability compared to purely data-driven approaches. One of the most widely used PIML frameworks is the physics-informed neural network (PINN) (Raissi et al., 2019), which augments standard neural network training with residual losses for differential equations, initial/boundary conditions, and observed data (Raissi et al., 2020). PINNs have demonstrated success across many fields, such as fluid mechanics (De Florio et al., 2021, 2022b; Cai et al., 2021), epidemiology (Schiassi et al., 2021a; Millevoi et al., 2024; Han et al., 2024; Kharazmi et al., 2021), cardiovascular mechanics (Sahli Costabal et al., 2020; Yin et al., 2021; De Florio et al., 2025), and beyond (De Florio et al., 2024; Ahmadi Daryakenari et al., 2024), while extensions incorporating variational inference, dropout, and deep ensembles have enabled principled quantification of aleatoric and epistemic uncertainty under noisy or incomplete observations (Hüllermeier and Waegeman, 2021; Zou and Karniadakis, 2025; De Florio et al., 2025).

Despite their flexibility, PINN‐based methods incur significant computational cost due to backpropagation through complex loss terms and often require careful tuning of weighting coefficients between physics and data objectives. To alleviate these challenges, the theory of functional connection (TFC) (Mortari, 2017b) provides an alternative “constrained expression” for differential equation solutions that analytically satisfy boundary or initial conditions (Mortari, 2017a). The resulting constrained approximation permits rapid, training‐free inference via direct least‐squares solves of residuals at collocation points (Schiassi et al., 2022). Building on TFC, the X-TFC (Schiassi et al., 2021b) employs a single‐layer random‐projection network whose input weights and biases are randomly sampled once and held fixed. This technique has been widely used recently with promising performance in terms of accuracy and computational efficiency in several PIML frameworks (Fabiani et al., 2021; Fabiani, 2025; Fabiani et al., 2024; Dong and Li, 2021; Wang and Dong, 2024; Dwivedi and Srinivasan, 2020; Galaris et al., 2021). X-TFC has been applied to both forward and inverse differential problems, demonstrating fast convergence and robustness in scenarios with limited, gappy, or noisy data (Osorio et al., 2025). Recent studies extend X-TFC to grey‐box and black-box system identification (De Florio et al., 2024; Ahmadi Daryakenari et al., 2024) and uncertainty quantification (UQ) in physiological models, showing its capacity to decompose total uncertainty into aleatoric, epistemic, and model‐form contributions (De Florio et al., 2025).

In this work, we review the conceptual foundations of PIML and X-TFC, apply X-TFC to an existing 1475 h vibration-based health-index (HI) dataset from a crack of the inner ring of a high-speed bearing in a 2.2 MW wind turbine gearbox (Bechhoefer and Dubé, 2020), and compare its effectiveness to classical physics-based fracture mechanics models. Finally, we show how our unified framework also successfully predicts the RUL of two bearings in a different gearbox model using an HI generated from supervisory control and data acquisition (SCADA) measurements (Eftekhari Milani et al., 2026), broadening its applicability to a wider range of industrial assets. Once the effectiveness of X-TFC has been assessed, future work will explore the application of X-TFC to other bearing cracking problems in wind turbines, such as pitch bearing ring cracks shown in Fig. 1c. We also employ Monte Carlo (MC) X-TFC to quantify epistemic uncertainty for different scenarios of dataset sizes. The paper is structured as follows: Sect. 2 describes the vibration-based HI dataset from the 2.2 MW wind turbine gearbox; Sect. 3 describes classical physics-based fracture mechanics models; Sect. 4 describes PIML and the X-TFC method; Sect. 5 assesses the effectiveness of X-TFC combined with a classical physics-based fracture mechanics model on the vibration-based HI dataset, using UQ analysis, and on two similar cracks with HIs derived from only SCADA data; and Sect. 6 summarizes the analysis.

2 Crack growth and health indicator dataset

As described in Sect. 1, the first dataset used in this paper is derived from an existing 1475 h (61 d) HI dataset from a crack of the inner ring of a high-speed bearing in a 2.2 MW wind turbine gearbox (Bechhoefer and Dubé, 2020). The HI represents a fusion of vibration-based condition indicators, collected at irregular intervals (approximately every 10 min on average, with occasional 20–30 min or multi-hour gaps) and used to extract features correlated with bearing damage. The HI for a normal bearing is approximately 0.1, while a value of 1 indicates that it is appropriate to do maintenance. The RUL is then the time from the current HI until the HI is 1. Note that this nomenclature does not define a probability of failure for the component, nor does it mean that the component fails when the HI is 1. Instead, it suggests performing maintenance before collateral damage or cascading faults occur. By performing maintenance on a bearing prior to shedding extensive material, costly gearbox replacement can be avoided, and the reliability of the gearbox can be restored to its design requirement (Bechhoefer and Dubé, 2020).

Figure 2a shows the HI calculated over the last 1475 h of operation of the wind turbine prior to the discovery of the bearing crack. Initially, the HI was very low; however, at 1000 h remaining, the HI began to rise, which corresponds to high loads from a winter storm (Bechhoefer and Dubé, 2020). By the time 400 h remaining was reached, the HI reached 0.5 on average, with maximum values greater than 1. At 0 h, the average HI exceeded 1, which triggered a borescope inspection of the bearing that uncovered the crack in the inner ring shown in Fig. 2b. It is not known whether this crack originated from WECs; however, it had the same axial orientation on the inner ring and occurred in a high-speed bearing like many WECs in wind turbine gearboxes. This HI dataset is used in the remainder of the paper to assess the effectiveness of X-TFC combined with a classical physics-based crack growth model.

Figure 2 Example vibration-based HI dataset and bearing crack. (a) Calculated vibration-based HI (Bechhoefer and Dubé, 2020). (b) Crack in 2.2 MW wind turbine gearbox bearing (Bechhoefer and Dubé, 2020).

3 Fracture mechanics models

As stated in Bechhoefer and Dubé (2020), fatigue crack propagation is commonly described in terms of three fundamental displacement modes (Beer et al., 1999; Frost et al., 1999): Mode I (opening), Mode II (in-plane sliding), and Mode III (out-of-plane shear). In Mode I, the crack faces separate under tensile loading; in Mode II, the faces slide laterally within the plane of the crack; and in Mode III, the faces move parallel to the crack front. By superposition of these modes, most general crack tip displacement fields can be represented. Under linear elastic assumptions, the crack surfaces are treated as traction-free boundaries so that external loads affect only the intensity of the stress field near the crack tip, which is quantified by the stress intensity factor K. The factor K depends on the component geometry and loading; for simple geometries, it scales linearly with the applied nominal stress σ and with the square root of the crack length a (Beer et al., 1999). For example, the Mode I stress intensity factor for a crack of length a can be expressed as

$$\begin{array}{}\text{(1)}& K=\mathit{\sigma }\sqrt{\mathit{\pi }a}\mathit{\alpha },\end{array}$$

where σ is the far-field (gross) stress, and α is a dimensionless shape factor. When K is known, stresses and strains near the crack tip can be calculated. For instance, the normal stress in the x direction at a point (r,θ) measured from the crack tip is (Frost et al., 1999)

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_x(r,\mathit{\theta })={\displaystyle \frac{K}{\sqrt{\mathrm{2}\mathit{\pi }r}}}\cos \left({\displaystyle \frac{\mathit{\theta }}{\mathrm{2}}}\right)\left[\mathrm{1}-\sin \left({\displaystyle \frac{\mathit{\theta }}{\mathrm{2}}}\right)\sin \left({\displaystyle \frac{\mathrm{3}\mathit{\theta }}{\mathrm{2}}}\right)\right],\end{array}$$

which illustrates the classical inverse-square-root singularity of the stress field. This formulation of the crack tip field in terms of a single parameter K is a cornerstone of linear elastic fracture mechanics. In fatigue crack growth analysis, the driving force is taken as the stress intensity range ΔK (the difference between maximum and minimum K in a load cycle). Experimental evidence (Frost et al., 1999) shows that ΔK, not the maximum K, governs crack advance, and if ΔK is held constant, then the crack growth rate becomes steady.

3.1 Linear elastic fracture mechanics

For many engineering materials (e.g., steels in gears and bearings) subjected to cyclic tensile loading, the incremental crack growth per cycle is Mode I, and it follows a power-law relation (i.e., Paris' law, Pugno et al., 2006). A commonly used form is

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}a}{\mathrm{d}N}}=D(\mathrm{\Delta }K{)}^m,\end{array}$$

where $\mathrm{d}a/\mathrm{d}N$ is the crack extension per load cycle, D is a material-specific constant, and m is the crack growth exponent (typically $\mathrm{3}\le m\le \mathrm{5}$). Substituting the expression $\mathrm{\Delta }K=\mathrm{2}\mathit{\sigma }\sqrt{\mathit{\pi }a}\mathit{\alpha }$ into Eq. (3) yields

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}a}{\mathrm{d}N}}=D{\left[\mathrm{2}\mathit{\sigma }\sqrt{\mathit{\pi }}\mathit{\alpha }\right]}^m{a}^{m/\mathrm{2}}.\end{array}$$

Rearranging Eq. (4) and integrating both sides give the total number of cycles N required to grow the crack from an initial length a₀ to some length a_f:

$$\begin{array}{}\text{(5)}& N=\underset{a_{\mathrm{0}}}{\overset{a_f}{\int }}{\displaystyle \frac{a^{-m/\mathrm{2}}}{D{\left(\mathrm{2}\mathit{\sigma }\sqrt{\mathit{\pi }}\mathit{\alpha }\right)}^m}}\mathrm{d}a.\end{array}$$

Performing the integration yields

$$\begin{array}{}\text{(6)}& N=-{\displaystyle \frac{\mathrm{d}N}{\mathrm{d}a}}{a}_{\mathrm{0}}-{\displaystyle \frac{a_f(a_{\mathrm{0}}/a_f{)}^{m/\mathrm{2}}}{m/\mathrm{2}-\mathrm{1}}}\end{array}$$

provided m≠2. In particular, if one assumes m=4, Eq. (6) reduces to

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}a}}=-{\displaystyle \frac{N}{a_{\mathrm{0}}\ln \left(\frac{a_f}{a_{\mathrm{0}}}\right)}}.\end{array}$$

Under steady, cyclic loading (e.g., constant-speed operations for helicopter gearboxes or for wind turbines in region 3), the number of cycles N is directly proportional to operating time. If the prognostic HI is assumed proportional to the accumulated damage (so that advancing crack length corresponds to decreasing HI), then Eq. (7) provides an RUL estimate by relating crack growth to the elapsed time or number of cycles. We can select a_f=1, defining the need to perform maintenance when the HI is 1.

3.2 Head's theory

Head's theory proposes a simplified model of the crack tip stress field, in which the material ahead of the crack tip is treated as an assembly of parallel elastic bars. Each bar of modulus E carries the applied remote stress σ both directly and through shear transfer. The resulting crack growth law can be cast in a form similar to Eq. (4) but with a higher effective exponent:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}a}{\mathrm{d}N}}={\displaystyle \frac{\mathrm{2}E{\mathit{\sigma }}^{\mathrm{3}}{a}^{\mathrm{3}/\mathrm{2}}}{\mathrm{3}{E}_W\left(\mathrm{\Delta }\mathit{\sigma }\right)D}}\propto K^{\mathrm{3}}.\end{array}$$

In effect, Head's model predicts a crack advance that corresponds to m=6 in the Paris law form. Although the intermediate steps involve material and geometric parameters, the final integrated form resembles Eq. (6) with m=6. In particular, inversion and integration analogous to the linear elastic case yield an alternate RUL expression:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}a}}={\displaystyle \frac{N}{a_{\mathrm{0}}\left(\mathrm{2}{a}_f-\frac{\mathrm{2}}{\sqrt{a_{\mathrm{0}}}}\right)}}.\end{array}$$

3.3 Dislocation theory

In the case of Mode III (anti-plane shear) dominant loading, dislocation (or plastic dislocation) theory may apply. Here the crack tip plastic zone is modeled as a continuous distribution of dislocations. Crack advance is assumed to occur when the integrated plastic slip reaches a threshold. The associated growth law can be expressed as

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{d}a}{\mathrm{d}N}}={\displaystyle \frac{a^{\mathrm{2}}{\mathit{\sigma }}_{\max }^{\mathrm{4}}}{DE{\mathit{\sigma }}^{\mathrm{3}}}}.\end{array}$$

This equation implies that crack advance per cycle grows with decreasing crack length, similar in form to the Paris law in Eq. (3), with an effective exponent m=2. By inverting, integrating, and changing terms, we get

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}a}}=-{\displaystyle \frac{N}{a_{\mathrm{0}}(\mathrm{2}{a}_f-\mathrm{2}\sqrt{a_{\mathrm{0}}})}}.\end{array}$$

3.4 Comparison of crack growth models

Figure 3 compares the three models described in the previous sections when applied to the full 1475 h vibration-based HI dataset described in Sect. 2. For this comparison, the values of a₀ and a_f are 0.05 and 1, respectively. With the given crack growth exponent m assumed in each model, Head's theory best matches the slow rise in the calculated HI data, especially over the last 700 h. In comparison, the linear elastic and dislocation theory models have a very slow rise, followed by a much steeper rise thereafter that would provide very little warning of the necessity of an impending maintenance action. Thus, in the remainder of this paper, we choose to use Head's theory model in combination with X-TFC.

Figure 3 Solutions of the differential equations given by linear, Head's theory, and dislocation theory models compared to the calculated vibration-based HI.

4 Physics-informed machine learning

Physics-informed machine learning (PIML) is an emerging framework that integrates known physical laws into the training of data-driven models. By embedding governing differential equations, boundary conditions, and other domain-specific constraints directly into the loss function, PIML enables models to learn from both observational data and prior physical knowledge. This approach contrasts with purely data-driven methods, which may struggle when observations are scarce, noisy, or incomplete.

Consider a physical system governed by partial differential equations,

$$\begin{array}{}\text{(12)}& \mathcal{N}\left[u\right(\mathit{x}),\mathit{x}]=\mathrm{0},\mathit{x}\in \mathrm{\Omega },\end{array}$$

where u(x) is the state variable, and 𝒩[⋅] is a differential operator defined on the domain Ω. Traditional data-driven models seek a surrogate $\widehat{u}(\mathit{x};\mathit{\theta })$ parameterized by θ that minimizes data loss:

$$\begin{array}{}\text{(13)}& {\mathcal{L}}_{\mathrm{data}}\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{d}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{d}}}\left|\widehat{u}\right({\mathit{x}}_i;\mathit{\theta })-u_i^{\mathrm{obs}}{|}^{\mathrm{2}},\end{array}$$

where $\mathit{\left\{}\right({\mathit{x}}_i,u_i^{\mathrm{obs}}){\mathit{\}}}_{i=\mathrm{1}}^{N_{\mathrm{d}}}$ values are N_d observed data points. In contrast, PIML augments this objective with a physics-based penalty:

$$\begin{array}{}\text{(14)}& {\mathcal{L}}_{\mathrm{physics}}\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}}}}\sum _{j=\mathrm{1}}^{N_{\mathrm{p}}}\left|\mathcal{N}\right[\widehat{u}({\mathit{x}}_j;\mathit{\theta }),{\mathit{x}}_j]{|}^{\mathrm{2}},\end{array}$$

where $\mathit{\{}{\mathit{x}}_j{\mathit{\}}}_{j=\mathrm{1}}^{N_{\mathrm{p}}}$ values are collocation points chosen in Ω. The physical loss is made by the residual of the differential equation, which needs to be minimized. In our case, it is the residual of Eq. (9):

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}a}}-{\displaystyle \frac{N}{a_{\mathrm{0}}\left(\mathrm{2}{a}_f-\frac{\mathrm{2}}{\sqrt{a_{\mathrm{0}}}}\right)}}=\mathrm{0}.\end{array}$$

The total loss is then given by

$$\begin{array}{}\text{(16)}& \mathcal{L}\left(\mathit{\theta }\right)={\mathit{\lambda }}_{\mathrm{data}}\cdot {\mathcal{L}}_{\mathrm{data}}\left(\mathit{\theta }\right)+{\mathit{\lambda }}_{\mathrm{physics}}\cdot {\mathcal{L}}_{\mathrm{physics}}\left(\mathit{\theta }\right),\end{array}$$

with weighting coefficients λ_data and λ_physics balancing data fidelity against physical consistency.

By enforcing Eq. (12) during training, PIML offers several advantages when dealing with realistic datasets:

• Regularization via physics. The physics-based term ℒ_physics acts as a principled regularizer, preventing overfitting caused by noise in $u_i^{\mathrm{obs}}$.

• Data efficiency. Incorporating prior knowledge reduces the reliance on large volumes of labeled data, enabling accurate learning even when N_d is small or spatial coverage is incomplete.

• Robustness to missing observations. The model can be evaluated at arbitrary collocation points, filling gaps in the observational domain and providing continuous field reconstructions.

• UQ. Extensions of PIML integrate Bayesian or ensemble approaches to quantify epistemic and aleatoric uncertainties, further improving reliability under data scarcity.

4.1 Extreme theory of functional connections (X-TFC)

The PIML framework used in this work is based on the X-TFC algorithm (Schiassi et al., 2021b), which differs from classic PINN for two main introduced features: the TFC (Mortari, 2017b, a; Mortari et al., 2019), which satisfies initial and boundary conditions of the differential equations analytically, and random-projection neural networks (Fabiani, 2025) trained via an extreme learning machine algorithm (Huang et al., 2006), thus avoiding the need for the computationally expensive backpropagation algorithm. The TFC constructs an explicit approximation of the solution of a differential equation that identically satisfies given boundary or initial conditions. For a general initial value problem of the form

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}=f\left(x\right(t),t),x\left(\mathrm{0}\right)=x_{\mathrm{0}},\end{array}$$

TFC represents the solution of the differential equation via a so-called “constrained expression”:

$$\begin{array}{}\text{(18)}& x(t;\mathit{\beta })=g(t;\mathit{\beta })-g(\mathrm{0};\mathit{\beta })+x_{\mathrm{0}},\end{array}$$

where g(t;β) is any sufficiently smooth function parameterized by coefficients β. In the X-TFC framework (De Florio et al., 2022a), g(t;β) is chosen as a random-projection single-layer neural network whose input weights and biases $\mathit{\{}{w}_j,b_j{\mathit{\}}}_{j=\mathrm{1}}^L$ are randomly sampled once and held fixed. Specifically,

$$\begin{array}{}\text{(19)}& g(t;\mathit{\beta })=\sum _{j=\mathrm{1}}^L{\mathit{\beta }}_j\mathit{\sigma }(w_j,t+b_j)={\left[\begin{array}{c}\mathit{\sigma }(w_{\mathrm{1}}t+b_{\mathrm{1}})\\ \mathit{\sigma }(w_{\mathrm{2}}t+b_{\mathrm{2}})\\ \mathrm{⋮}\\ \mathit{\sigma }(w_{L}t+b_L)\end{array}\right]}^T,\end{array}$$

and

$$\begin{array}{}\text{(20)}& {\mathit{\sigma }}^T\left(\mathrm{0}\right)=\left[\begin{array}{c}\mathit{\sigma }\left(b_{\mathrm{1}}\right)\\ \mathit{\sigma }\left(b_{\mathrm{2}}\right)\\ \mathrm{⋮}\\ \mathit{\sigma }\left(b_L\right)\end{array}\right]={\mathit{\sigma }}_{\mathrm{0}}^T,\end{array}$$

where σ(…) is the neural network's activation function, and β is the vector of the output weights of the neural network to be learned.

Thus, a parametric approximation of the solution of the ordinary differential equation and its derivative in time are

$$\begin{array}{}\text{(21)}& x(t;\mathit{\beta })={\left[\mathit{\sigma }-{\mathit{\sigma }}_{\mathrm{0}}\right]}^T\mathit{\beta }+x_{\mathrm{0}}\end{array}$$

and

$$\begin{array}{}\text{(22)}& \dot{x}(t;\mathit{\beta })=\dot{\mathit{\sigma }}(t{)}^T\mathit{\beta }.\end{array}$$

A least-squares formulation minimizes the residuals of the ordinary differential equation at the collocation points {t_i} and any available data. Denoting the vector of ordinary differential equation residual loss

$$\begin{array}{}\text{(23)}& {\mathcal{L}}_{\mathrm{physics}}={\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}-f\left(x\right(t),t)\end{array}$$

and data loss

$$\begin{array}{}\text{(24)}& {\mathcal{L}}_{\mathrm{data}}=x(t;\mathit{\beta })-\widehat{x}(t;\mathit{\beta })\end{array}$$

(where $\widehat{x}(t;\mathit{\beta })$ represents the observation data) by $\mathcal{L}={\mathcal{L}}_{\mathrm{physics}}+{\mathcal{L}}_{\mathrm{data}}$, the update step in a Gauss–Newton scheme to determine the kth iterate of the unknown coefficients β is

$$\begin{array}{}\text{(25)}& \mathrm{\Delta }{\mathit{\beta }}^k=-[{\mathcal{J}}^T\mathcal{J}{]}^{-\mathrm{1}}{\mathcal{J}}^T\mathcal{L},\end{array}$$

where 𝒥 is the Jacobian of losses with respect to the coefficient vector β. A schematic of the X-TFC framework is shown in Fig. 4. X-TFC offers rapid, training-free inference by exploiting random projections and a direct least-squares solver, making it well suited for online estimation and UQ.

Figure 4 Schematic of the X-TFC algorithm that includes data and physics in the loss function. Input weights and biases are randomly selected. The last step solves the least-squares problem.

5 Results

As previously mentioned, the physics loss consists of the residual of the differential equation modeled by Head's theory as a fracture mechanics model, where the values of a₀ and a_f are 0.05 and 1, respectively. The simulations are run with 1000 collocation points of the neural network's input, sampled equally spaced between −1 and 1, and with five hidden neurons. The numbers of collocation points and neurons are user-selected to achieve a good compromise between accuracy and computational expense. The input weights and bias are sampled with a uniform distribution between −1 and 1, and the selected activation function is tanh. Figure 5 reports the results for all our simulations. Each plot represents a different data size scenario for different RUL estimates to understand the advantage of using PIML, especially when data availability is scarce, such as when the crack is just beginning to form, or sparse, such as when the HI is less frequently calculated. For each data size scenario, we performed the RUL prediction for different weights given to the differential equation, such as ${\mathit{\lambda }}_{\mathrm{physics}}=[\mathrm{0.0},\mathrm{0.25},\mathrm{0.50},\mathrm{0.75},\mathrm{1.0}]$, where scenario λ_physics=0 is the case in which we do not have knowledge of physics, so we only perform data-driven regression. All simulations were performed with λ_data=1.

The scatter data points in blue represent the full dataset we have for the bearing under study, while the red scatter data points are those used in our framework, simulating conditions of scarce data and early stage crack growth. The time domain, which represents the RUL in hours, is calibrated to correspond to zero when the HI reaches a value of 1, which means that there is a high probability of bearing damage and a need to perform maintenance, obtained with the scenario of full knowledge of physics (λ_physics=1) and full availability of data (see Fig. 5a). This gives us our ground truth solution. The results in Fig. 5a represent the case in which we have all data throughout the entire crack growth period, allowing us to see the effect of the weights on the physics loss to the RUL prediction. When we assume no knowledge of physics, the data alone are sufficient to have a good prediction of RUL, which differs from our ground truth by about 56 h, as reported in Table 1.

Table 1 Error values of RUL prediction in hours for different combinations of data availability and physics weights. In the right column, the computational time in milliseconds is reported.

Figure 5 Results of PIML for different dataset sizes. For each scenario, different weights (${\mathit{\lambda }}_{\mathrm{physics}}=[\mathrm{0.0},\mathrm{0.25},\mathrm{0.50},\mathrm{0.75},\mathrm{1.0}]$) on the equation loss are given to understand the improvement in the RUL forecasting brought by the physics knowledge. (a) Full dataset available: 100 %. (b) Partial dataset available: 75 %. (c) Partial dataset available: 50 %. (d) Partial dataset available: 25 %.

Figure 6 Error in RUL prediction (h) for different scenarios of dataset size and weights on physics loss.

The plots in Fig. 5b, c, and d show a linear trend of the deviation of the RUL prediction, as the availability of data and knowledge of physics decrease, as expected. This is clear evidence of the improvement brought by physics that reduces the search space previously based only on data. In fact, it can be observed that even with low data availability (25 %), the RUL prediction has an error of only about 54 h when full weight (physics weight = 1.0) is given to the loss of the differential equation. However, with no contributions from physics (physics weight = 0.0), the RUL predication error increases to ≈600 h. In Fig. 6, the RUL error surface is plotted, obtained with 9191 simulations (91 dataset sizes from 0.1 % to 100 % and 101 physics weights from 0 to 1), showing a mirror-like behavior that explains the equal importance of both data and physics, fundamental ingredients of PIML, and the advantage of using both of them in synergy to obtain the maximum performance gain. With data availability even as low as 10 %, a physics weight greater than 0.8 keeps the RUL error under 60 h. The simulations in this work were performed on a MacBook Pro with an Apple M3 Pro (12 cores, 36 GB) using MATLAB R2024b, resulting in the computational times for each case reported in Table 1.

5.1 Uncertainty quantification

To quantify epistemic uncertainty, we employed a simple yet scalable MC strategy, which is called MC X-TFC (De Florio et al., 2025). In this method, we independently train multiple X-TFC models on synthetic data corrupted by random noise drawn from a known distribution, and each model starts from a different random initialization of weights and biases. MC X-TFC closely parallels the “deep ensemble” frameworks, as proposed in Lakshminarayanan et al. (2017), which has proven to be highly effective for neural network UQ, even when that uncertainty arises solely from random parameter initialization (Ganaie et al., 2022; Rahaman and Theiry, 2021; Wenzel et al., 2020; Gawlikowski et al., 2023; Psaros et al., 2023; Zou et al., 2024b; Pickering et al., 2022; Zou et al., 2024a; Zhang et al., 2024). We define epistemic uncertainty as the variability in the predicted state vector that arises from two sources:

1. model‐parameter uncertainty, that is, the fluctuations in the coefficients β induced by the random sampling of network weights w and biases b in Eq. (19)

2. data‐driven uncertainty, that is, the variability in the observed data under our noise model, which we also propagate into the initial condition a₀ of the HI by perturbing it with zero‐mean Gaussian noise with standard deviations of 20 %.

In our analysis, we evaluate the predictive performance of the model using the mean of error (ME), the standard deviation of error (SDE), and the confidence interval (CI) of signed errors, where errors retain their direction (positive or negative). All error values and derived statistics are expressed in hours, and the simulations are performed with 100 MC realizations and physical loss weight equal to 1. The ME is calculated as

$$\begin{array}{}\text{(26)}& \mathrm{ME}={\displaystyle \frac{\mathrm{1}}{m_c}}\sum _{i=\mathrm{1}}^{m_c}{e}_i,\end{array}$$

where e_i (in hours) is the error for the ith sample (predicted minus true value), and m_c is the total number of samples. The ME indicates systematic bias: a positive ME means earlier predictions of damage and the need for maintenance, and a negative ME means later predictions. The SDE is given by

$$\begin{array}{}\text{(27)}& \mathrm{SDE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{m_c-\mathrm{1}}}\sum _{i=\mathrm{1}}^{m_c}(e_i-\mathrm{ME}{)}^{\mathrm{2}}},\end{array}$$

which measures the spread of signed errors around the ME. A larger SDE (in hours) indicates greater variability in the prediction deviations. We also compute non-parametric percentiles of the e_i distribution to capture coverage analogous to Gaussian confidence levels:

• 68 % interval (“1σ”), 16th to 84th percentiles of e_i;

• 95 % interval (“2σ”), 2.5th to 97.5th percentiles of e_i;

• 99.7 % interval (“3σ”), 0.15th to 99.85th percentiles of e_i.

These intervals indicate where the central 68 %, 95 %, and 99.7 % of signed errors lie, respectively, preserving information on whether prediction errors tend to be positive or negative. All these error metrics for epistemic UQ are reported in Table 2 for different dataset sizes. As expected, the absolute values of ME and SDE increase with a decrease in dataset size. The negative values of ME represent an error in an overly early prediction of damage and the need for maintenance. Computational times also decrease with a decrease in dataset size because of the smaller number of points in the data loss function. The CIs reported in Table 2 are also qualitatively shown in Fig. 7 for the three confidence percentages. Despite the inherent stochasticity of random model initialization and noisy observations, the resulting CIs remain moderate (e.g., ±10 h for a 68 % CI on the full dataset). This behavior stems primarily from the strong physics constraint given by Head's theory differential equation, which tightly regularizes the solution and limits divergence across ensemble members. As such, once model-induced variability is largely suppressed by physics-informed regularization, the spread in predictions stays relatively small, highlighting the effectiveness of embedding physical laws into the learning process.

Figure 7 Epistemic quantification for different dataset sizes, considering λ_physics=1 on the equation loss, and three CIs with 67 %, 95 %, and 99.7 % confidence obtained with 100 MC realizations. (a) Full dataset available: 100 %. (b) Partial dataset available: 75 %. (c) Partial dataset available: 50 %. (d) Partial dataset available: 25 %.

Table 2 Error values of RUL prediction in hours for different combinations of data availability while quantifying epistemic uncertainty. The error metrics in hours are the mean of error (ME), the standard deviation of error (SDE), and the confidence interval (CI) of error. The last column reports the computational time in seconds for 100 MC realizations.

Figure 8 Ablation study on epistemic quantification for different random initializations of input weights and bias ranges, considering λ_physics=1 on the equation loss; partial data size (25 %); and three CIs with 67 %, 95 %, and 99.7 % confidence obtained with 100 MC realizations. Random initialization of input weights and bias in ranges of (a) [−1, 1], (b) [−10, 10], and (c) [−20, 20].

Figure 9 Aleatoric uncertainty band estimated from the vibration-based HI. The shaded region represents ±2σ bounds based on the time-varying standard deviation of the residuals between noisy data and smoothed data. This illustrates the heteroscedastic, irreducible uncertainty arising from measurement noise.

Since the epistemic uncertainty is reducible, the CIs presented in Fig. 7 are obtained with optimized parameters that reduce the epistemic uncertainty as low as possible. This reduction depends on the many hyperparameters selected, such as the activation function, number of collocation points, number of neurons, least-squares iteration tolerance, and random sampling ranges. An ablation study on the behavior of epistemic uncertainty depending on the choice of the random initialization of the neural network input weights and bias is shown in Fig. 8. The three plots show that the different choices of input weights and bias ranges can lead to wider CIs for the RUL prediction, and by choosing the optimal hyperparameter setup, the epistemic (model) uncertainty can be reduced. Further reduction in prediction spread would therefore require enhancements in data quality, such as using higher-resolution measurements, to minimize noise-driven variability. Moreover, the aleatoric uncertainty is inherently irreducible (Hüllermeier and Waegeman, 2021) and, in this case, is particularly high due to the very noisy data. In practical terms, this means that there is a minimum level of predictive uncertainty, driven by data noise, which cannot be further reduced by improving the model. Consequently, even a perfect model would still exhibit this baseline uncertainty dictated by the natural variability of the observations. This is shown in Fig. 9, where the CI of the aleatoric uncertainty is plotted against the noisy dataset. Because the dataset exhibits heteroscedastic noise, that is, noise amplitude varies over time, we first estimate the global noise scale by analyzing the full HI series. We apply a locally estimated scatter plot smoothing (LOESS) method to find a smoothed trend in the HI over time and then compute the residuals using the noisy dataset to capture fluctuations. Then, the overall noise level is approximated as the standard deviation of the residuals. Finally, the aleatoric uncertainty band in Fig. 9 is generated by plotting the smoothed trend ±2σ, reflecting the heteroscedastic, irreducible noise-driven uncertainty in the HI measurements.

Figure 10 Epistemic uncertainty quantification with SCADA-based HI of two different gearbox bearings, considering λ_physics=1 on the equation loss, and three CIs with 67 %, 95 %, and 99.7 % confidence obtained with 100 MC realizations. (a) Gearbox bearing no. 1. (b) Gearbox bearing no. 2.

5.2 Application to SCADA-based health index

To demonstrate the generality and adaptability of the X-TFC framework, we further apply our methodology to an HI derived from SCADA data described in Eftekhari Milani et al. (2026). The SCADA dataset is from 1.5 MW wind turbines, with measurements available in a 10 min time frame beginning at their commissioning. These turbines have the same gearbox model, with one experiencing an intermediate-speed shaft bearing axial crack and the other experiencing a high-speed shaft bearing axial crack. The cracks are the first damage that occurred to these gearboxes, enabling an analysis of the bearing condition from a pristine state to cracked.

In the data pre-processing step, the non-operational data points with zero rotor speed or negative power are removed. In addition, the time stamps with missing or erroneous values due to sensor problems are dropped. Then, a daily averaging is performed, decreasing the missing data percentage from 14 % in the 10 min time frame to around 2 % in the daily time frame. This averaging increases data continuity and reduces the computational burden. To obtain the HIs, the approach proposed in Eftekhari Milani et al. (2026) is used. First, for each turbine, the gearbox bearing temperature in the healthy state is modeled by training a convolutional neural network using the remaining six SCADA channels during the first year of operation, in which the bearings are expected to be healthy. Then, the modeled healthy-state gearbox bearing temperature values are subtracted from the recorded values, obtaining the residuals. They are then normalized according to the 1-year training period to obtain the residual z scores. These z scores, along with the four operational and environmental signals of power, rotor speed, wind speed, and ambient temperature, are input to a convolutional autoencoder (CAE)-based unsupervised HI construction model to build the HIs. A leave-one-out cross-validation is performed, where one wind turbine is set aside to test the CAE, and the remaining undamaged wind turbines are used to train the CAE. The trained CAE is then used to obtain the test wind turbine's HI.

While purely data-driven methods like LSTM and XGBoost showed promise for SCADA data, their precision was limited, often resulting in a high number of false alarms (Desai et al., 2020). To address this limitation, we applied X-TFC to the SCADA-based HIs of the two turbines with cracked bearings (Eftekhari Milani et al., 2026), embedding the same Head's theory differential equation used previously for vibration-based data. The physics model acts as a regularizer, enforcing consistency between the SCADA-derived HI trajectory and the expected fatigue crack growth behavior. In Fig. 10, epistemic uncertainty in the X-TFC RUL prediction for the SCADA-data-based HI is reported. When trained on only 25 % of the available SCADA HI trajectory, X-TFC obtains failure predictions with a confidence level of about 10 d. Importantly, the physics-informed structure of X-TFC helped to suppress model overfitting compared to traditional machine learning models applied directly to SCADA signals (Desai et al., 2020). This confirms the versatility and applicability of the X-TFC framework to both vibration- and SCADA-data-based HIs. This widens the scope of prognostic applications for wind turbines, enabling early and reliable failure prediction even when only typical SCADA data are accessible.

6 Conclusions

In this work, we have introduced and demonstrated the application of X-TFC as a PIML framework for predicting the RUL of wind turbine gearbox bearing cracks. By embedding Head's fracture mechanics model directly into the learning objective, X-TFC balances fidelity to HIs based on both vibration and SCADA data with consistency to the differential crack growth law. In particular, we calibrated the baseline PIML RUL estimate from the full 1475 h vibration-based HI dataset from a cracked high-speed bearing in a 2.2 MW wind turbine gearbox and full physics weighting (λ_physics=1) to a value of HI equal to 1, indicating a high possibility of bearing damage and a need to perform maintenance. Lower values of weights and purely data-driven regression (λ_physics=0) still lead to good accuracy in RUL prediction, with less than an error of 60 h. Even when only 25 % (or less) of the vibration-based HI data are available, heavy weighting on the physics loss yields RUL errors below 60 h, a four-fold improvement over the data-only case (≈600 h). This demonstrates X-TFC's ability to leverage domain knowledge to fill gaps in noisy, sparse, or intermittent monitoring data, while requiring very low computational times (of the order of milliseconds). To further enhance the interpretability and robustness of these predictions, we introduced an MC ensembling strategy (MC X-TFC) to quantify epistemic uncertainty. This approach produces CIs for the estimated RUL by propagating uncertainty due to both model initialization and measurement noise, thereby enabling probabilistic forecasting. These intervals provide insight into the reliability of each prediction and help operators assess risk in real-time maintenance planning. Finally, as shown for this case, the irreducible aleatoric uncertainty is high due to the natural variability of the HI data.

These results demonstrate the potential of X-TFC as a powerful, data-efficient, and computationally light methodology for bearing fatigue crack prognostics in rotating machinery. By guaranteeing physically plausible degradation trajectories regardless of data sparsity, X-TFC enables earlier and more reliable maintenance decisions, with the potential to reduce unplanned downtime and O&M costs for wind turbine gearboxes. Our study, in its early stage, focuses on HI data derived from bearing cracks and employs Head's theory as the governing crack growth law. By validating our methodology using HI data generated from both vibration and SCADA datasets, we confirm its versatility to bearings of multiple gearbox manufacturers, models, and ratings. To capture more realistic wind turbine bearing behavior, our PIML framework should be validated against larger, more diverse datasets. These models would enable real-time monitoring of bearing health by continuously ingesting streaming data and updating predictions as new information arrives.

Future extensions should also examine fatigue cracking in pitch bearing rings, which may require more comprehensive fracture mechanics models, such as multi-phase Paris laws or dislocation theory blends, to capture more complex and realistic nonlinear crack–load interactions. Continued development of UQ capabilities, such as MC X-TFC, will also be critical for producing probabilistic RUL estimates with actionable confidence bounds, especially in high-stakes decision environments.