Enabling Safe Fast Charging of Lithium-Ion Batteries via a Simulation-Trained Digital Twin

1. Introduction

With the rapid and large-scale adoption of electric vehicles (EVs), increasing battery energy density and reducing charging time have become primary development objectives. Among these, fast charging has gained significant attention as an effective means to shorten charging duration. However, high-rate charging is well known to accelerate battery degradation, leading to reduced capacity and increased internal resistance [1,2]. Especially in aged cells, fast charging generates substantial heat which can be difficult to dissipate effectively, thereby exacerbating degradation mechanisms and thermal stress within the cell [3], which may bring cell operation into a dangerous zone where battery thermal runaway (TR) triggering is more likely [4]. Li et al. [5] comparatively analyzed the effects of different charging rates on thermal runaway characteristics and demonstrated that increasing the charging rate significantly reduces battery thermal stability, which highlights the strong coupling between fast-charging conditions, thermal behavior, and safety risks in lithium-ion batteries.

Predicting failure in lithium-ion cells is therefore crucial for ensuring the safety of EVs and other battery-powered applications. Accurate early prediction of thermal runaway would enable preventive measures to be implemented before critical conditions are reached. When mitigation strategies are activated in a timely manner, thermal runaway propagation can be prevented, thereby minimizing damage to the battery system and the vehicle [6] and potentially to people. Therefore, defining a safe operating area in which batteries function with a minimal likelihood of thermal runaway is a critical task.

Fault prognostics in battery management systems (BMSs) typically rely on model-based or data-driven approaches. Model-based methods [7,8,9,10], employ electrochemical, electrical, thermal, or coupled physical models to estimate battery states; however, their scalability and computational cost represent a challenge, especially for large battery systems. Furthermore, ageing effects during operation are hardly taken into consideration. In contrast, data-driven approaches utilize machine learning and statistical techniques to analyze operational data—such as voltage, current, and temperature—along with environmental conditions and historical performance to detect faults and predict failures [11,12,13,14,15,16,17]. These methods depend on experimental measurements, simulation data [16,18], or publicly available datasets, such as NASA Ames [19], CALCE [20], and the Stanford/MIT fast-charging datasets [21], which have been utilized across a wide range of applications [22,23,24,25,26,27,28,29,30,31,32] including state-of-health prediction, fault diagnostics and remaining useful life (RUL) estimation.

A digital twin is a real-time virtual representation of a physical object or system, used to monitor, analyze, and simulate its behavior and performance. Traditional Battery Management Systems rely on simple rules and often ignore complex electrochemical dynamics, which can lead to reduced accuracy and reliability, especially under demanding conditions like high temperatures or fast charging [33]. Therefore, coupling digital twin with BMS can enhance accuracy and reliability. Recent studies have explored digital twin frameworks for lithium-ion batteries, primarily focusing on state estimation [34,35], performance monitoring [33], and degradation prediction [36,37]. However, despite their advancement, these approaches are largely deterministic and do not explicitly address rare but safety-critical events such as thermal runaway, nor do they incorporate probabilistic risk quantification under uncertain operating conditions.

On the other hand, most studies focus on defining the Safe Operating Area (SOA) for battery operation, which typically outlines permissible voltage, temperature, and current ranges to avoid degradation [38,39]. Research has demonstrated that fast charge cycling can reduce thermal stability and lower the thermal runaway onset threshold, while high charging rates combined with battery ageing increase thermal hazard conditions that conventional SOA limits may not adequately address [40]. Recent work on early warning models also suggested that safety assessments based on temperature alone may be insufficient, underscoring the need for more comprehensive guidelines and predictive tools to ensure safe battery operation in practical fast-charging scenarios [41].

To address these gaps, in this work a simulation-trained digital twin framework is proposed that integrates high fidelity electrochemical–thermal modeling, machine learning, and probabilistic risk analysis to enable the prediction of thermal runaway conditions and the estimation of critical charging current during fast charging, effectively connecting fault prognosis with safety-oriented decision support. By explicitly linking operating conditions to TR proneness-based critical charging currents, the proposed approach enables risk-aware fast-charging limits identification that would balance performance and safety. The key contributions of present study are:

(i) the use of an encoder–decoder architecture enabling efficient surrogate modeling of static-to-dynamic Li-ion cell thermal behavior;

(ii) the formulation of thermal runaway proneness as a probabilistic metric using Monte Carlo sampling under parameter uncertainty; and

(iii) the introduction of a risk-based critical charging current, which links acceptable failure probability thresholds to actionable operating limits.

2. Methods

The methodology adopted for this study consisted of 3 parts: 1) creation of a cell digital twin for temperature calculations and thermal runaway occurrence prediction, 2) estimation of quantitative thermal runaway proneness for specific cell condition, and 3) calculation of the relative critical charging current. In the following sections, the methodology adopted for each part of the study, will be discussed in detail, while in Section 3 an example of application will be presented.

2.1. Digital Twin

2.1.1. Digital Twin Creation

As mentioned above, a simulation-trained digital twin is defined as a model representing thermal behaviour of a Li-ion cell. For data generation, a battery model was set up in COMSOL Multiphysics v6.3 environment [42], based on the Single Particle (SP) electrochemical-thermal coupled approach, which assumes the electrodes as spherical particles and a uniform current distribution across the battery [43]. This model had been validated against data on NCA and NMC 18650 cylindrical cells in an ambient with natural convection for predicting temperature and thermal runaway behavior of the cells [7].

To simplify the analysis and allow a more direct comparison, constant cell chemistry, ambient temperature and heat transfer conditions were preliminarily imposed, i.e. NCA/graphite electrodes, 24 °C initial temperature and natural convection conditions on the cell surface. It is worth noting that different chemical compositions and environmental parameters can be studied, along with more complex electrochemical models, adopting the same framework established with the current analysis.

The influence of three fundamental parameters has been here investigated, i.e. charging current with CC regime, cell capacity and internal resistance, by developing a script, written in Python using MPh open-source library v1.3.1 which uses COMSOL’s native Java APIs to automate simulation runs. The results of the simulations consisted of cell state of charge and temperature temporal profiles. In order to allow a direct comparison, each run was continued until the cell reached a full (100%) state of charge.

For all simulations the model calculated the temperature increase with respect to ambient temperature over time, within the cell. In Figure 1 the cell average temperature deviation is reported for two sample cases, both conditions triggered thermal runaway.

Within data-driven framework, since the input features are static while the output (temperature trend) is dynamic, conventional Artificial Neural Networks (ANN) or standalone Long Short-Term Memory (LSTM) architectures are not well suited to directly model this relationship. In fact, a simple ANN treats the output sequence as a fixed vector and does not explicitly capture temporal dependencies, while a standard LSTM is designed for sequence-to-sequence mappings and is not naturally suited for static-to-dynamic relationships. Differently, an encoder–decoder architecture addresses this limitation by encoding the static inputs into a latent representation, which is then transformed into a sequence and processed by a decoder part into a time-dependent temperature trajectory. Synthetically, an encoder–decoder architecture is a neural network design where one part, the encoder, reads and compresses input data into a meaningful internal representation, while the second part, the decoder, handles that representation and generates a final dynamic output [44,45].

Accordingly, an encoder–decoder architecture was here adopted to effectively capture the temporal dynamics of the cell thermal behaviour, and its schematic visualization is presented in Figure 2. In particular, the encoder consisted of a feedforward artificial neural network with two hidden layers, which map the static input features of the system into a latent representation layer. The latent representation is then passed through a so-called RepeatVector layer, which assumes the same information at each time step and replicates it across the 50 time-steps in which the whole charge phase (i.e. state-of-charge, SOC, increasing from 0 to 100%) is subdivided. This transforms the static encoded features into a sequence suitable for the decoder. LSTM networks are then effectively employed in the decoder stage. The decoder processes the latent sequence to generate the temperature evolution over time.

Despite the representational power gained from deeper configurations, such models are susceptible to overfitting, which occurs when a model accurately learns the specific patterns of the training data without understanding the basic underlying relationships, thus leading to poor performance in unseen data representation. To mitigate that effect, a dropout rate of 25% of units was applied on each layer.

The set-up model was trained and evaluated using the temperature profiles generated from COMSOL simulations, with the whole dataset randomly split into 80% for training and 20% for validation.

Training machine learning models requires selecting appropriate “hyperparameters” (external variables that configure the training process), which need to be tuned to optimize model performance. Hyperparameters tuning involves selecting optimal training, typically by minimizing metrics such as a loss function. In this study, hyperparameters were tuned using Grid Search, a brute-force approach that evaluates all combinations within a predefined set, which is suited for small-scale problems. The hyperparameter search space included the number of units (32 and 64), dropout rate (0.1 and 0.25), and activation function (ReLU and tanh). The choice of activation function governs the nonlinear transformation applied at each neuron. ReLU (Rectified Linear Unit) passes positive pre-activations unchanged while setting negative values to zero, introducing sparsity and accelerating convergence, and the hyperbolic tangent function (tanh) squashes outputs to the range [−1, 1], providing a zero-centred alternative that can be beneficial in recurrent layers.

2.1.2. Model Evaluation

The mean absolute error (MAE) and the root mean squared error (RMSE) were chosen here to evaluate the model performance against the simulations-generated dataset; MAE measures the magnitude of errors without considering their direction, making it less sensitive to outliers, while RMSE quantifies the magnitude of prediction errors giving greater weight to larger deviations, thus making it sensitive to outliers and calculation inaccuracies.

$$\mathrm{M}\mathrm{A}\mathrm{E}={\displaystyle \frac{1}{\mathrm{n}}}\sum \left|{\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}}\right|,$$

(1)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}\sum {({\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}})}^2},$$

(2)

In the above equations, ${\mathrm{y}}_{\mathrm{i}}$and $\widehat{{\mathrm{y}}_{\mathrm{i}}}$ are the values obtained from the electrochemical-thermal simulations and those predicted by the DT model, respectively, and n is the number of runs. Model development and evaluation were carried out in Python open-source API of Keras and Scikit-learn.

Overall, the Digital Twin predicted the temperature profiles with reasonable accuracy, with an MAE of 3.05 °C, and an RMSE of 6.69 °C. Figure 3 illustrates the encoder-decoder model results and the COMSOL simulation outputs for three different runs. As it can be seen, the encoder-decoder model performed rather well in predicting the overall temperature trend and, most importantly, it predicted with a high level of reliability thermal runaway occurrence under the different conditions investigated.

After developing and evaluating the performance of the DT model, it is necessary to assess its reliability outside the training domain. This step is essential to assess whether the temperature profiles generated by the DT remain physically meaningful when applied to unseen conditions. Therefore, an analysis was conducted to further explore the physical interpretability of the DT outputs, under a set of selected conditions, all at a current rate of 2C (Table 1).

2.1.3. Thermal Runaway Assessment

Beside predicting the temperature trend within the cell, the main aim of this study was to assess the possibility of thermal runaway occurrence under given operating conditions, and this can be accomplished once the onset condition of the sequence of side reactions leading to TR has been preliminarily defined. With reference to this, it is commonly recognized that the first side reaction preceding the thermal runaway sequence is the decomposition of the Solid Electrolyte Interphase (SEI) [46,47,48], which provides the heat required for initiating the subsequent exothermic reactions. Consequently, the onset temperature of the SEI decomposition was adopted as a threshold to indicate the occurrence of thermal runaway. As reported in the literature [46], this temperature depends on multiple conditions such as cell chemistry and electrolyte composition. Since in this study the digital twin was trained on an NCA/graphite cell chemistry, the onset temperature was assumed to be about 80 - 90 °C [46,47,48,49], corresponding to an increase of around 60 °C from the initial ambient temperature (see also Figure 1).

Development of a digital twin requires significant preliminary modeling effort and computational resources, therefore, its usage for thermal runaway prediction compared to other simpler data-driven approaches must be properly motivated. An alternative Artificial Neural Network (ANN) was here adopted to evaluate and compare the capacity of a purely data-driven model in TR prediction. This comparison enabled assessing whether the digital twin provided added value in terms of accuracy, robustness, and generalization capability. The ANN shared the DT feature encoder architecture, but its output layer consisted of a single neuron to reflect the binary classification objective, namely whether a thermal runaway would occur under the given conditions.

To ensure a consistent basis for comparison, the ANN was trained on the same dataset used for the DT model, and the results of the comparison between the two approaches will be presented in the study case, based on their relative confusion matrix and on the following related parameters:

• Accuracy: the percentage of all predictions that were correctly classified.
• Precision: the ratio of true positive to all positive predictions made.
• Recall: the proportion of the actual positive cases that were correctly identified.

2.2. Cell Proneness to Thermal Runaway

Thermal runaway of a cell is very unlikely when each of the operating parameters lies well within its safe range. Conversely, the closer a parameter is to its safe operating limits, the higher the likelihood of incurring in a thermal runaway; in such a condition, even small variations in its value, for example due to control issues, sensitivity level of sensors or instrumentation accuracy, can bring the cell working into a hazardous area. Prediction of these conditions would therefore help preventing cell operation in such a dangerous “pre-failure” area. As an attempt to solve this issue, a procedure has been devised to calculate what might be defined the “proneness” of the cell to thermal runaway, in a given condition characterized by the nominal values of the operating parameters.

However, due to the complex and interacting effects of the operating parameters, each of them characterized by specific safe limits and differently influencing cell response, determining a safe area is not a straightforward task. In order to get a more realistic result, in the present study all the three investigated parameters (current intensity, cell capacity and resistance) were considered simultaneously. These three parameters define a three-dimensional parameter space, and each operating condition of the cell is represented by a point within this space. Considering that each sensor of a BMS is subject to unavoidable measurement uncertainty [50], each parameter is associated with a variability range. Accordingly, a continuous population centred around the nominal operating point can be defined. Within this volume, some conditions will lead to thermal runaway, while others will not. The proportion of conditions leading to thermal runaway over the entire population, is defined here as the “proneness” of the cell to thermal runaway under the considered nominal conditions. In order to estimate this parameter, a random sample was drawn from this population, and by using the developed DT model of the cell, the conditions that resulted in thermal runaway were identified. Subsequently, the proneness of the cell to thermal runaway was quantified by the frequency of occurrence of thermal runaway in the sampled volume as:

$${\mathrm{P}}_{\mathrm{T}\mathrm{R}}={\displaystyle \frac{{\mathrm{n}}_{\mathrm{T}\mathrm{R}}}{{\mathrm{n}}_{\mathrm{t}}}},$$

(3)

where P_TR, n_TR and n_t are proneness to thermal runaway, the number of conditions resulting in a thermal runaway and the total number of runs investigated (sample size), respectively.

In order to ensure an unbiased and systematic exploration of the parameters space, a Monte Carlo simulation method was introduced, with the input parameters randomly sampled assuming a normal distribution for each of them, around their specific nominal values. Since it is recognized that the reliability of any statistical approach is influenced by the sample size, an appropriate value must be determined to minimize its effect. To determine an appropriate sample size, the proneness to thermal runaway was calculated for a range of sample sizes with fixed input conditions and ranges. The minimum sample size at which the proneness converged, i.e. it showed negligible variations with further increases in sample size, was finally selected.

In order to understand the influence of standard deviations (variability range) of the studied parameters on proneness calculation, a 4-level full factorial design was then used adopting the values reported in Table 2.

Finally, the dependency of the proneness to TR on the input parameters were also investigated. Therefore, a multiple linear regression analysis was carried out, adopting normalized independent variables (current, capacity and resistance). The resulting coefficients provided a measure of the influence of each factor on the proneness to thermal runaway. This sensitivity analysis was performed using a full factorial design, with the values reported in Table 3 and using Python with open-source API of pyDOE and sklearn libraries.

2.3. Critical Current Calculation

By setting an arbitrary reference value for proneness to thermal runaway, it would be possible to identify operating regions that should not be exceeded to avoid entering conditions where the likelihood of experiencing a thermal runaway would become unacceptable. This might be done with reference to any of the operating parameters and in particular, from a fast charging perspective, to the charging current, with CC regime.

To calculate the critical current, a threshold for thermal runaway proneness must first be defined, above which operation of the cell is considered unsafe: in this study, four reference thresholds, i.e. 0.01, 0.05, 0.1 and 0.2 were adopted.

To determine the critical current, the problem is formulated as a root-finding task. Specifically, a function was defined as the difference between the predicted thermal runaway proneness and a predefined safety threshold, and the critical current would correspond to the value for which this function equals zero, as described in Equation 4:

$$\mathrm{O}\mathrm{b}\mathrm{j}=\mathrm{P}\left(\mathrm{I}|\mathrm{Q},\mathrm{R}\right)-{\mathrm{P}}_{\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}}$$

(4)

where Obj, P and P_threshold are the objective function, the proneness as a function of current at specific capacity and resistance, and the proneness threshold previously specified. Brent’s method [51] is then employed to locate the root of the objective function. This method combines bisection (halving technique) and interpolation techniques, providing both robustness and efficiency. It is particularly well suited to this application because it guarantees convergence when a valid bracket is available while requiring relatively few function evaluations, which is important given the computational cost of evaluating the DT-based proneness model.

3. Study Case

The methodology described in Section 2 is here demonstrated on a commercially representative 18650 NCA/graphite cell, simulated under constant-current charging at an ambient temperature of 24 °C. The study case is structured in three stages. First (Section 3.1), the physical interpretability of the digital twin was assessed by examining its cell thermal behaviour prediction capability across a systematic set of conditions outside the training domain, as anticipated in Section 2.1.2. Second, the Monte Carlo proneness framework is applied to quantify thermal runaway risk as a function of operating parameters and their associated uncertainties (Section 3.2). Finally, the critical charging current is calculated and tracked over the ageing life of battery no. 5 from the NASA dataset [19], demonstrating how the framework can inform risk-aware fast-charging limits in a realistic degradation scenario (Section 3.3).

3.1. Thermal Behavior

To investigate the physical meaning of the DT model, Figure 4 presents the simulated thermal behavior of a NCA 18650 cell under several operating conditions (see Table 1) using the digital twin model, developed according to the described methodology. Figure 4a explores the effect of varying cell resistance (respectively, 10, 30, and 60 mΩ) at fixed capacity (Q=0.5 Ah) and current (I=2C). The temperature deviation from the initial ambient temperature, clearly increases with resistance increase, with the highest resistance (60 mΩ) leading to a substantially sharper rise and a higher final temperature compared to the lower resistances. It must be observed that the duration of the charge phase reported in Figure 4 is not in absolute terms but in time steps, as described in the methodology section, to get an immediate comparison of all the simulation runs under the different conditions. As expected, Figure 4a indicates that higher internal resistance significantly contributes to heat generation and temperature rise in the battery during operation.

Figure 4b shows the temperature deviation for fixed resistance R=30 mΩ and current I=2C, while varying the capacity Q from 0.5 Ah to 0.8 Ah, and 1 Ah. The plot demonstrates that as capacity increases, the temperature deviation rises more rapidly and reaches higher maximum values. The results shown in this figure are consistent with previously reported studies [7].

These results suggest that the digital twin model is able to provide physically meaningful simulations of the thermal behaviour of cells also under conditions different than those adopted for training of the model itself.

In order to check the capability of the DT model in predicting TR occurrence, a confusion matrix is reported in Figure 5a, comparing DT predictions with the COMSOL simulations; in the matrix, 1 and 0 indicate a positive or negative TR occurrence, respectively. It can be seen that the DT was capable of identifying the majority of thermal runaway events, missing only 2 out of 72 occurrences, i.e. with a recall of 97.22 %; similarly, it properly identified 115 cases where TR did not occur, out of 120 total runs, with a ratio of more than 95.8 %. The overall accuracy (percentage of correct predictions) reached 96.35 %.

As described in the methodology section, these values were then compared in Table 4 with the equivalent ones obtained with an alternative artificial neural network model adopted as a benchmark (see Figure 5b). In the latter case, a recall and an accuracy of 93 % and 94.8 % were obtained, respectively, denoting an average improvement for the recall of more than 4% with the digital twin approach. This result, along with the possibility of generating full temperature details, in contrast with the ANN approach, support the use of the DT model for further analyses.

3.2. Proneness to Thermal Runaway Analysis

The devised digital twin has been adopted for calculating cell proneness to thermal runaway (PTR) under a number of operating conditions. As discussed in Section 2.2, this requires a preliminary definition of the study volume, which depends on some system characteristics, such as sensors sensitivity, etc., which will define the range of variability (standard deviation) of each study-parameter. In the present study case, standard deviations of 0.05, 0.2 Ah and 10 mΩ for C-rate, capacity and resistance, were used, respectively.

Figure 6 presents the dependency of PTR on the sample size (total number of simulations performed) for three cell operating conditions. The plot shows how the proneness evolves as the number of simulations increase. It can be seen that for small sample sizes, the estimate is highly unstable, denoting a condition of insufficient sampling. As the number of runs increases beyond roughly 1,000, the estimate begins to settle down, and from about 10,000 onward only slighter fluctuations are observed. Those remaining fluctuations are small and appear to be random rather than systematic, indicating that the estimator has effectively converged. Choosing 50,000 runs as a sample size, seems a safe and reasonable compromise because it lies well within the stable region of the curve, far beyond the initial transient phase where larger changes occur. Further increasing the number of runs would involve an unnecessary much higher computational effort without significant beneficial effect.

Figure 7 presents the sensitivity of proneness to thermal runaway to the standard deviations of all three input parameters, organised as three pairs of panels — each pair fixing one parameter variability at a high and a low level while varying the remaining two — for a nominal operating condition of I = 1C, Q = 3 Ah and R = 70 mΩ.

In Figure 7a,d, the current standard deviation was fixed at its high (0.06C) and low (0.015C) values, respectively. In Figure 7a, proneness increases consistently with R_std across all Q_std levels. The highest capacity variability (Q_std = 0.15 Ah) produces a markedly steeper and higher curve — reaching PTR ≈ 0.018 at R_std = 5 mΩ — while the three lower Q_std values cluster closely together in the range 0.010–0.012, indicating that the additional contribution of capacity variability becomes significant only at its highest level. In Figure 7d, by contrast, it can be observed that reducing current variability to 0.015C suppresses proneness by roughly an order of magnitude across the entire parameter space. Only the highest capacity variability (Q_std = 0.15 Ah) produces a noticeable rise at large R_std, reaching a maximum proneness of just about 0.00065; all lower Q_std curves remain essentially flat near zero. This comparison establishes that when current variability is small, resistance and capacity variability contribute negligibly to thermal runaway risk.

In Figure 7b,e, capacity standard deviation was fixed at 0.15 Ah and 0.011 Ah, respectively. In Figure 7b, the four I_std curves are separated: I_std = 0.06C reaches P_TR ≈ 0.018, I_std = 0.045C reaches P_TR=0.0075, I_std = 0.03C remains at 0.0025, and I_std = 0.015C is indistinguishable from zero across the full range of resistance variability. In Figure 7e, where capacity variability is reduced to its minimum, the same qualitative ordering is preserved but absolute proneness values are lower, with P_TR at I_std = 0.06C reaching ≈ 0.010 and at I_std = 0.045C reaching about 0.003, with reductions of approximately 45% and 60% respectively relative to the high-capacity-variability case. This confirms that capacity variability provides a meaningful secondary contribution, but only when current variability is already at a moderate-to-high level.

Figure 7c,f represent cases with fix resistance standard deviation at 5 mΩ and 1 mΩ, respectively. In Figure 7c, proneness rises with Q_std for all I_std levels, with I_std = 0.06C beginning at PTR ≈ 0.010 even at the lowest capacity variability and reaching about 0.018 at Q_std = 0.15 Ah, while lower I_std curves rise more gradually. In Figure7f, halving the resistance variability substantially reduces proneness across the range — the maximum proneness (I_std = 0.06C, Q_std = 0.15 Ah) falls to ≈ 0.009 — and, importantly, for I_std ≤ 0.03C proneness remains near zero regardless of how large capacity variability becomes. This indicates that the contribution of capacity variability to proneness is not independent: it requires a minimum threshold of current uncertainty to become detectable.

Overall, Figure 7 consistently identify current variability as the primary driver of thermal runaway proneness. Resistance and capacity variability act as secondary amplifiers whose contributions scale with I_std and become meaningful only when current uncertainty is sufficiently large. This hierarchy is reflected directly in the standardized regression coefficients reported in Table 5, deriving from the sensitivity analysis carried out on the three input parameters and already introduced in Section 2.2 (Table 3).

The standardized coefficients reported in Table 5 provide a direct measure of the relative influence of each input parameter on the response, independently of their original units. Among the considered factors, current exhibited the highest standardized coefficient value, indicating that it has the strongest influence among the analysed parameters on the likelihood of thermal runaway. Capacity follows with a coefficient of 0.4202, suggesting a moderate contribution, while internal resistance showed a slightly lower coefficient value, indicating a comparatively smaller, yet still meaningful, effect. While all three parameters contribute to thermal runaway proneness, current appears to be the dominant factor within the investigated parameter space. This result has an important practical implication: reducing current uncertainty or variability (e.g. by improving the accuracy of current sensing in a BMS) has the greatest reduction effect on thermally induced risk, particularly under fast-charging conditions where operation lies closer to the thermal runaway boundary.

3.3. Critical Current Estimation

The concept of proneness to thermal runaway has finally been applied to the calculation of the critical current during charge, once a threshold value of the proneness had been set for the cell under investigation. In order to check the methodology, battery no. 5 in the NASA dataset [19] was used, a 18650 commercial Li-ion battery, for which detailed information during the charging/discharging tests, all run at ambient temperature of 24 °C, was available.

Since the method in this study required knowledge of the internal resistance as an input parameter, while the NASA dataset did not report this specific value, a parameter estimation had to be preliminarily performed. Figure 8a illustrates the capacity trend over cycles number reported in the NASA dataset, while Figure 8b presents the corresponding resistance, calculated from simulations [52]; the progressive degradation behaviour of the cell with ageing is apparent from both figures.

Determining the critical (maximum) current depends on the proneness threshold adopted, so that, in its essence, the plot quantifies the trade-off between safety (low thermal runaway proneness) and performance (high allowable current) for Li-ion cells. Different thresholds were here assumed, and the results are illustrated in Figure 9: as the proneness threshold increases (which translates into a higher tolerance for thermal runaway risk) the critical current also increases, clearly indicating that the cell is allowed to operate at higher currents when more risk is accepted. The three curves, related to randomly selected cycles in NASA dataset, show similar trends demonstrating that, although ageing or cycle variations slightly shift the absolute critical current, the overall relationship remains consistent.

Another result which can be derived from Figure 9 is that the maximum C-rate trend cannot be easily predicted, but specific calculations are required to properly estimate it. In fact, looking at the three curves, it can be observed that, at equal risk tolerance, the maximum C-rate initially increases with cycling ageing (compare cycles 52 and 125), but then it decreases (from cycle 125 to 149).

These considerations can also be seen in Figure 10, where the evolution of the critical current (at 0.10 proneness threshold) for battery no. 5 from NASA dataset is reported. As it can be seen, the critical current starts increasing with ageing, primarily driven by capacity fade, but, after an almost stationary phase, it decreases due to the effect of increased resistance. In fact, although a continuous decrease would be expected with ageing based on the continuous rise in cell internal resistance, capacity fade remains the dominant factor governing thermal runaway susceptibility in the early stages [7,53], as also represented in Table 5. As degradation progresses and the cell approaches the capacity knee point, the influence of resistance becomes increasingly significant, leading to a shift in the governing mechanism of critical current behavior [53].

The possibility of characterizing the critical current is a useful result also from a practical point of view, since it would allow identifying the maximum charge current (and corresponding minimum charge time), taking into account the true history of the cell (cycling ageing), i.e. its actual parameters values over time, and not only based on its nominal properties. At the same time, risk tolerance can also be included as a parameter in operation conditions choice.

4. Conclusions

This study presents a simulation-trained digital twin framework that delivers a risk-calibrated critical charging current for fast-charging operation, defined as the maximum allowable charging rate that maintains the probability of cell thermal runaway below an application-defined threshold. Unlike conventional detection-based approaches, which identify hazardous conditions after onset, this framework enables predictive determination of safe operating limits that can be directly integrated into CC fast-charging strategies. In order to do this, a proneness metric is introduced to quantify the deviation of operating conditions from the hazardous operating area, providing a systematic measure of thermal runaway risk.

An encoder–decoder neural network was developed to map static operating conditions (charging current, capacity, and internal resistance) to dynamic temperature evolution. Assuming a temperature-based criterion to characterize the onset of thermal runaway sequence, a probabilistic framework is introduced using Monte Carlo simulation to quantify cell proneness to thermal runaway under parameter uncertainty. Among the investigates parameters charging current was identified as the dominant factor, followed by capacity and internal resistance, confirming that current control is the most impactful lever for fast-charging safety management.

Since the model was trained entirely on simulation data its accuracy is bounded by the adopted assumptions so that effects such as localized thermal gradients and cell-to-cell variability are not explicitly captured. Also, the uncertainty ranges assigned to the input parameters were representative values for framework demonstration rather than derived from specific experimental characterization.

Nonetheless, it is believed that the proposed approach can represent a useful tool for accident prevention and operational optimization. Future work should focus on validating the framework against experimental thermal runaway data across multiple cell chemistries, incorporating multi-factor uncertainty distributions informed by real BMS sensor accuracy, and extending the approach to multi-cell pack configurations where thermal interactions between cells introduce additional risk pathways. Integration of the critical current estimator within an online BMS loop, where ageing state is continuously updated, represents the most direct pathway toward practical deployment of this framework in fast-charging safety management.