Concept of Digital Twins (DTs)

The concept of digital twins (DTs) has gained significant traction in recent years, particularly in the engineering and industrial disciplines [1]. A DT refers to a virtual model that closely mirrors a real-world physical product, system, or process, serving as an effectively indistinguishable digital counterpart for practical purposes such as simulation, integration, testing, monitoring, and maintenance. DTs have been recognized as the fundamental premise of product lifecycle management, encompassing the entire lifecycle of the physical entity they represent, including creation, construction, operation/support, and disposal [2-6].

The value and fidelity of a DT representation depend on the specific use cases, considering that its granularity and complexity are tailored accordingly. Notably, a DT can be conceptualized and utilized even before its physical counterpart exists, allowing for comprehensive modeling and simulation of the intended entity’s lifecycle.

The origin of the DT concept can be traced back to NASA, which pioneered its practical definition in 2010 as part of its efforts to enhance the simulation of spacecraft using physical model representations [2, 6]. Since then, the development of DTs has progressed in tandem with advances in product design and engineering. Traditional manual drafting techniques have given way to computer-aided drafting and design and model-based system engineering approaches, which establish a direct link between the DT and signals from its physical counterparts.

In recent years, DTs have garnered significant attention in various industrial sectors. Notably, their application in the chemical processing industry has demonstrated versatility [7]. In addition, DTs are increasingly utilized in simulation-based vehicle certification and fleet management [8], as well as in addressing the challenges of unpredictable and undesirable behaviors in complex systems [9, 10]. A specific area of growth has been observed in the nuclear plant industry, where digital twin architectures have been developed for the management and monitoring of nuclear plants, reflecting the expanding scope and potential of technology [11, 12].

Previous Work on Reactor Operation Digital Twins (RODTs)

A significant amount of research has been conducted to explore the potential applications of reactor operation digital twins (RODT) in the nuclear energy domain. RODTs represent a specific instantiation of the DT concept, focusing on the numerical representation of nuclear reactors for real-time prediction, optimization, monitoring, control, and decision-making during their operational stage [13].

A previous study [13] introduced the first prototype of an RODT specifically tailored for the prediction of neutron flux and power fields in the operational stage of the HPR1000 reactor core [14]. The prototype demonstrated the feasibility of using the RODT for online monitoring and prediction of key reactor parameters, thereby enhancing operational understanding and performance assessment.

In [15], the forward model of the RODT was realized through the application of a nonintrusive reduced-order model constructed using an SVD autoencoder to learn the nonlinearity of the field distribution. Machine learning techniques, specifically leveraging k-nearest neighbour (KNN) and decision trees, are employed to establish forward mapping. Additionally, the inverse model adopts a generalized latent assimilation method for accurate estimation of the model parameters.

While demonstrating promise, opportunities remain for improving the efficacy of RODT for practical deployment. In particular, the inverse solver relies on a computationally expensive methodology that limits its potential integration within an online framework. Therefore, the aim of this work [16] is to enhance the key capabilities of the RODT. This includes proposing an advanced differential evolution algorithm to upgrade the inverse solver for improved efficiency and accuracy of parameter estimation. In addition, uncertainty quantification was conducted to validate the RODT’s performance considering noisy observational data. Numerical validation experiments were performed across representative domains for an HPR1000 pressurized water reactor core to demonstrate its potential for engineering applications.

The key steps for constructing an RODT include: (i) training a reduced-order model (ROM) using model order reduction (MOR) [17-20] techniques as well as an efficient forward model and (ii) adaptation and implementation of inverse models to infer the input parameters and the resulting field distribution.

Introduction to Reduced Order Models (ROMs)

Many mathematical models encountered in real-life processes present computational challenges when employed in numerical simulations because of their inherent complexity and large dimensionality. To address these challenges, MOR techniques have been developed to reduce the computational complexity of these problems, particularly in the context of simulations involving large-scale dynamic and control systems. The ROM was computed as an approximation of the original model by reducing the dimensionality or degrees of freedom associated with the model.

The preparation of an ROM model involves two distinct phases, the offline and online phases, which are integral to the creation of a simplified representation of a high-fidelity full-order numerical model, enabling efficient real-time simulations and analysis. The offline phase focuses on determining the inherent low-dimensional structure of the underlying full-order model and deriving the reduced basis functions from high-fidelity snapshots of the system. Intrusive MOR techniques [21] are employed when a comprehensive understanding of the governing equations and numerical strategies employed in the full-order model is available. This requires access to the detailed numerical framework of the full-order model, which may involve proprietary commercial codes. By projecting a full-order model onto a reduced space, intrusive ROMs provide a concise representation with reduced computational complexity.

However, in practical engineering scenarios [22-24], detailed information regarding a full-order model is often unattainable, and the solvers or codes implementing these models are treated as black-box entities. This limitation poses a challenge to the implementation of traditional intrusive MOR methods. To address this issue, nonintrusive MOR techniques have been developed. Non-intrusive MOR approaches aim to construct ROMs without relying on detailed knowledge of the numerical framework of a full-order model. Instead, these methods establish an input-output mapping between the input parameters and the reduced basis through data-driven techniques such as interpolation, regression, or machine learning, approximating the behavior of the full-order model based on a reduced set of parameters and basis functions.

ROMs prepared through the offline and online phases play a critical role in the development of RODTs. RODTs require the integration of high-fidelity models and solvers with real-time capabilities, and ROM-based emulators provide promising solutions [21-23, 25]. Leveraging ROM-based forward solvers, RODTs can perform both one-to-one forward simulations and one-to-many simulations for inverse problems, encompassing parameter identification [26-29], data assimilation [24, 30-35], sensitivity/uncertainty quantification [36-40], and optimization/control [37, 41-47].

Forward Model Construction

In line with the framework RODT, as investigated in a previous study [13, 15], we selected KNN [48] to construct the forward model. The choice of KNN is justified by its numerous advantages, which render it a suitable candidate for the forward model in this context.

First, it is a nonparametric algorithm, meaning that it does not assume any specific functional form for the relationship between the inputs and outputs. Second, KNN is relatively simple to implement and computes efficiently, particularly for low-dimensional input problems. The complexity of the algorithm depends primarily on the number of training samples, making it suitable for problems with moderate dataset sizes [49]. In a recent study [50], an investigation was conducted to explore the possibility of utilizing different values of K for each dataset, considering the information provided by the correlation matrix.

Inverse Model Construction

Generally, the inverse problem can be divided into two categories:

• Parameter Identification. Solving the parameter vector based on the information from observation vectors: In this scenario, we aim to estimate the unknown parameters of the system based on the available observation vectors. The observation vectors contain information about the system’s behavior or response under conditions that are characterized by the parameter vector.

• State Estimation. Solving the field distribution based on the information from observation vectors: In this scenario, the objective is to determine the unknown field distribution of the nuclear system based on the available observation vectors.

Contribution of this Work

A previous work [13, 15, 16] demonstrated the feasibility of employing the RODT framework in engineering applications. However, a need exists for a systematic and in-depth investigation of the core components of RODT, namely parameter identification and state estimation. In addition, the development of mature and modular software is crucial to facilitating the application of RODT in practical scenarios. This motivated us to develop an efficient and appropriate surrogate inverse model to accurately predict the parameters and states. This also motivated our research on hybrid optimization approaches that can effectively balance coarse- and fine-grid search (coarse–fine-grid search) methods, which are discussed in the following sections.

The contributions of this research are as follows:

• Investigation and implementation of novel gradient-free optimization algorithms specifically tailored to handle the challenges posed by the discontinuous surrogate forward mapping constructed with the KNN forward model of our RODT and comparing the convergence profile, stability, and accuracy-time performance.

• Exploration and evaluation of the proposed coarse–fine-grid search methods within various global optimization approaches, aiming to strike a balance between accuracy and computational efficiency in parameter identification and state estimation.

• Development and deployment of efficient modular RODT software integrated with the aforementioned algorithms into the inverse problem solver, which applies to other sustainable energy systems.

These contributions focus on advancing state-of-the-art RODT by addressing the unique challenges encountered in solving inverse problems in nuclear engineering and proposing novel methodologies for improved parameter identification and state estimation. Furthermore, our work addresses the critical challenge of achieving real-time, high-precision simulations in nuclear energy DT applications, significantly enhancing the capabilities of DTs in the nuclear energy sector and facilitating improved monitoring, control, and optimization of nuclear energy systems.

Offline phase: Training a Non-Intrusive ROM

The variation in the power field Φ in a nuclear core is characterized by physical laws that are generally described implicitly by a governing equation [51], such as the neutron transport or diffusion equations, which is written as Equation 1. $\{\begin{array}{l}\mathcal{F}\left(\Phi \left(r\mathrm{,}\mu \right)\mathrm{,}\mu \right)=0\hfill \\ r\in \Omega \subset {ℝ}^d\hfill \\ \mu \in \mathcal{D}\subset {ℝ}^p\hfill \end{array},$ (2.1) where $\text{Ω}\subset {ℝ}^d$ represents the spatial domain of dimension d, with $d\ge 1$, and $\text{ϕ}\left(\mu \right)=\text{Φ}\left(r\mathrm{,}\mu \right)$ is defined in a Hilbert space that is equipped with an inner product $\langle \cdot \mathrm{,}\cdot \rangle$, and the induced norm $\Vert \cdot \Vert =\sqrt{\langle \cdot \mathrm{,}\cdot \rangle }$. $\mathcal{D}\subset {ℝ}^p$ represents the p-dimensional feasible parameter domain that covers the operation of the reactor core. In this study, Eq. 1 represents two-group neutron diffusion equations [52].

2.1.4

Training a Forward Model

The goal of the forward model is to construct an efficient surrogate mapping model between the parameter and coefficient, as follows: $\alpha \left(\mu \right)\cong F^{\text{ML}}\left(\mu \right)\mathrm{,}$ (2.9) where F^ML denotes the surrogate forward mapping constructed by the machine learning approach, such as model-dependent machine learning methods or model-free optimization searching methods.

The response of the ROM model, denoted as $F^{\text{ML}}\left(\mu \right)$, corresponds to the predicted coefficient data. The image of 𝑫 under the forward map $F^{\text{ML}}\left(\mu \right)$ represents the set of responses for all possible states. The difference $Y_{\text{o}}-\mathcal{H}Q{F}^{\text{ML}}\left(\mu \right)$ is referred to as the observation data misfit or the residuals associated with parameter μ.

We design the cost function as Eq. 12 based on the residual. The forward model is used in the search of the optimal approximation for the operator ϕ in the functional set $\left\{Q{F}^{\text{ML}}\mathrm{,}{F}^{\text{ML}}\in \mathcal{E}\left(D\mathrm{,}{ℝ}^n\right)\right\}$, where $\mathcal{E}\left(D\mathrm{,}{ℝ}^n\right)$ denotes the set of functions mapping from 𝑫 to ${ℝ}^n$, with n the reduced dimension in the MOR process.

Previous research [13] evaluated various forward model approaches, among which the KNN algorithm demonstrated exceptional accuracy and efficiency. In this study, we specifically selected the KNN algorithm [48] as the forward model to exploit its full potential. To enhance the performance of the KNN algorithm, we focused on two crucial parameters that significantly influence its efficacy.

The first parameter relates to the voting distance in KNN, wherein we investigate the impact of utilizing either the Euclidean distance or Manhattan distance as the metric for determining the nearest parameter choice in the training set. Notably, these two distances can be expressed as Minkowski norms with an exponent of p = 1, 2, corresponding to the Euclidean and Manhattan distances, respectively. The Minkowski distance, a measure of the discrepancy between the testing and training parameters, can be mathematically represented as follows: ${\Vert {\mu }^{\text{test}}-{\mu }^{\text{train}}\Vert }_p={\left({\displaystyle \sum _{i=1}^{dim\left(\mu \right)}{\left|{\mu }_i^{\text{test}}-{\mu }_i^{\text{train}}\right|}^p}\right)}^{\frac{1}{p}}\mathrm{,}$ (2.10) where $dim\left(\mu \right)$ denotes the dimension of the parameter μ. Finally, we choose the Euclidean metric.

Another vital parameter is K, which refers to the number of candidates voted for. The optimal choice of K can be determined by plotting the relative L₂ error of the predicted observation and true vectors, which is expressed as follows.

From Fig. 2, we can deduce that when K = 5, our forward model using KNN performs with the best prediction accuracy, and the variation in the error is relatively smooth, which corresponds to the continuous variation in the value of the field at the corresponding position of the field manifold snapshot.

Fig. 2Full-screen View

Finding the optimal choice of K for KNN model (Forward on the left, Inverse on the right).(a) KNN Forward model; (b) KNN Inverse model

In contrast, the stair-shaped observation data predictions in Fig. 3 demonstrate that the KNN forward surrogate model has a step-like progression as neighbors are adjusted between discrete datasets and thus are not continuous or differentiable. In this manner, well-developed state-of-the-art algorithms based on gradient information for solving optimization and inverse problems did not work well in our case. To meet this demand, we investigated other novel approaches, which are discussed in detail in the following sections.

Fig. 3Full-screen View

Non-differentiable stair-shaped observation data prediction via the KNN surrogate forward model (tested when the first parameter St varies from 0 to 615 steps, the other parameters take the midpoint values of the parameter intervals, respectively, 1250 for Bu; 50 for Pw; and 295 for Tin). (a) L₂ norm of predicted observation; (b) 1th vertical level, the 43th assembly; (c) 18th vertical level, the 11th assembly

Online phase: Parameter Identification and State Estimation

The inverse problem in nuclear engineering refers to the task of determining the model parameters that yield the observed data, as opposed to the forward problem, which involves predicting the data based on the given model parameters. The objective is to determine the model parameters μ^* that approximately satisfy Equation 12.

The nature of the inverse problem depends on whether operator ϕ is linear or nonlinear. In most cases, particularly when solving systems governed by the neutron flux diffusion function, the inverse problem is nonlinear, owing to the nonlinearity of the forward map. Thus, we attempt to develop nonlinear algorithms, gradient-free algorithms, or algorithms based on the initial guess given by fast linear algorithms, such as KNN.

Exhaustive Direct Search with Latin Hypercube Sampling (LHS)

To evaluate and compare the different approaches to solving inverse models, it is essential to establish a benchmark for measuring their performance. In this section, we propose using exhaustive direct search (EDS) with latin hypercube sampling (LHS) as the benchmark method for inverse problems. It combines the accuracy of EDS methods, which have been widely used since the 1960s [60], with the sampling efficiency of LHS.

LHS [61] is a sampling technique that ensures a comprehensive exploration of the solution space while maintaining desirable properties. The range of each variable was divided into equally probable intervals, and the sample points were strategically placed to satisfy the requirements of a Latin hypercube, where each sample was the only one in its corresponding axis-aligned hyperplane. The advantage of LHS lies in its independence from the number of dimensions, which allows for efficient sampling without the need for an increasing number of samples. Furthermore, LHS facilitates sequential sampling, enabling the tracking of the already selected samples.

Exhaustive Direct Search with 2 Steps Latin Hypercube Sampling (LHS2STEPS)

To enhance the benchmark algorithm and mitigate the uncertainties introduced by the randomness of the sampling process, [13] introduced a two-step sampling strategy. This strategy involves further sampling the neighborhoods of the n₁ best candidates generated during the initial stages of the LHS algorithm. For the detailed pseudocode, please refer to 3.3, excluding the part related to ALHS.

Exhaustive Direct Search with Assembled 2 Steps Latin Hypercube Sampling (ALHS)

To improve our benchmark algorithm and further reduce the uncertainty caused by the randomness of the sampling process, we introduce a mean mechanism at the end of the classical LHS algorithm to calculate the mean of the n₁ top-rated candidates in the five clones, which are spread around the n₂ top-rated candidates in the first sampling process. The concept is illustrated in Algorithm 1.

Efficient Deterministic Machine Learning Algorithm: KNN

In addition to the benchmark LHS method, another efficient deterministic algorithm that can be considered for solving inverse models is the KNN algorithm. KNN is a nonparametric classification and regression algorithm that can be adapted to solve inverse problems.

The KNN algorithm operates based on the principle of similarity. Given a new input data point, KNN finds the K-nearest neighbors in the training dataset based on the Euclidean distance and then generates the output value by averaging the output values of the K-nearest neighbors. The selection of the optimal value of K can be determined by graphing the relative L₂ error between the predicted observation vector and the true vector, as shown in Fig. 2. The analysis in Fig. 2 reveals that selecting K = 1 yields the optimal outcome in terms of minimizing the L₂ prediction error, thereby emphasizing its advantageous accuracy.

Metaheuristics Algorithms with Physical Constraints

Recently, the application of metaheuristic algorithms with physical constraints [62-64] has gained significant attention for inverse modeling. These algorithms offer powerful gradient-free optimization techniques that can effectively handle complex inverse problems while satisfying the physical constraints imposed by system dynamics.

In this manner, we implemented and adapted practical algorithms and constructed related solvers that were integrated into the developed RODT framework. These solvers are designed to address the inverse problem presented by Equation 12.

Hybrids of Optimization Methods for Inverse Problem: Balancing Coarse–Fine-Grid Search

To address the challenges posed by the discontinuity of KNN surrogate forward mapping and noisy observations, we propose hybrid algorithms that combine the strengths of deterministic machine-learning methods for coarse-grid search and metaheuristic algorithms for fine-grid search. This methodology can also be found in fields that encompass feature extraction of aerial data [77], automatic free-form optics design [78], and hyperparameter optimization [79]. The objective of a coarse–fine-grid search is to strike a balance between global exploration and local exploitation [80] to achieve accurate, efficient, and less sensitive parameter estimation. The hybrid optimization hub designed accordingly can be further illustrated in Fig. 4.

Fig. 4Full-screen View

Hybrid optimization schema of coarse–fine-grid search

3.6.2

Advanced CSDE (ACSDE) and Advanced DE (ADE)

We utilize the SVC algorithm to classify the test parameters and denote the center of the classified hypercube as the initial machine learning estimate μ_SVC; then, we use the CSDE and DE algorithms to further explore the parameter space with the initial population spread around the center μ_SVC. In the SVC process, according to a previous research [15] on the importance rate of the input parameters, the last parameter, temperature Tin, contributes slightly, and we also consider that Tin varies at a relatively small boundary; thus, we only administer the SVC process to the first three parameters. We chose two representative metrics to evaluate the performance of the SVC model.

Confusion Matrix. The confusion matrix [82] provides a comprehensive overview of the binary classification results by displaying the number of true positive (TP), true negative (TN), false positive (FP), and false negative (FN) predictions. This allowed us to assess the accuracy of the model’s by differentiating classes and identifying potential sources of misclassification.

We can infer from Fig. 5 that our SVC classification model predicts the first two parameters with relatively high accuracy while exhibiting confusion concerning the third parameter.

Fig. 5Full-screen View

(Color online) Confusion Matrix Plot that evaluates the efficacy of SVC prediction in the first 3 parameters (St, Bu, Pw, respectively)

F1-score. The F1-score [83], defined in Eq. 3.18, is a popular performance metric for binary classifier models. In Eq. 3.18, TP, FN, FP are the numbers of true positives, false negatives, and false positives classified by the model, respectively. $F1=\frac{TP}{TP+\frac{1}{2}\left(FP+FN\right)}$ (3.18) The predictive model employed in this study for each parameter is a multi-class SVC utilizing the “one-vs-others" strategy. This approach involves training separate binary classifiers to distinguish between the candidate intervals for a given parameter. Consequently, performance evaluation necessitates the adoption of the Macro F1-score as a metric. The Macro F1-score (Eq. 19) is computed as the average of the F1-scores obtained for each interval, reflecting the overall effectiveness of the classifiers in capturing intra-parameter variations: $\text{F}{1}_{\text{macro}}:=\frac{1}{N}{\displaystyle \sum _{i=0}^N\text{F}{1}_i}\mathrm{,}$ (3.19) where i is the interval index, and N is the number of intervals.

Based on the Confusion Matrix shown in Fig. 5 and $\text{F}{1}_{\text{macro}}$, we designed the grid size of the fine-grid search process as follows:

Setup of Reactor Physics Operation Problems

Setup of background information for HPR1000. The objective of constructing the RODT was to accurately predict the power distribution within the HPR1000 reactor core during operation. The core consists of 177 vertical nuclear fuel assemblies, 44 of which are equipped with self-powered neutron detectors (SPNDs) to measure neutronic activity and power fields. Figure 6 provides a visualization of a horizontal slice of the HPR1000 core and an axial slice of an SPND-equipped assembly, with only one-quarter of the core displayed owing to the symmetry along the x and y axes. The gray fuel assemblies represent those containing control rods, whereas the assemblies marked with D indicate the presence of SPNDs. For more detailed information regarding the HPR1000 reactor and generic neutronic physical model, please refer to [14]. Further descriptions of the model for data assimilation and the initial work on RODT can be found in [35] and [13], respectively.

Fig. 6Full-screen View

Quarter of the core in the radial direction (white square: fuel assembly with SPNDs, grey square: fuel assembly with control rods, D: fuel assembly with neutron detectors) [13]

During the normal operation of the HPR1000 reactor, two types of control rods were utilized for regulation. The first type, known as compensation rods, is responsible for coarse control and a substantial reactivity reduction. These compensation rods comprise four subtypes: G1, G2, N1, and N2. The second type, called regulating rods (R rods), are employed for fine adjustments to maintain the desired power or temperature [84]. Power evolution within the reactor is influenced by various factors, including the movement of the control rods, the burnup of nuclear fuel, variations in the power level of the reactor core, and fluctuations in the inlet coolant temperature. This evolution was mathematically modeled using two-group diffusion equations [13] and numerically solved using the CORCA-3D code package [85]. Developed at NPIC, the CORCA-3D code is capable of solving 3D few-group diffusion equations, considering thermal-hydraulic feedback, and performing pin-by-pin power reconstruction.

Setup of parameter domain. In this study, we consider CORCA-3D to be an opaque solver in which we input the parameter vector μ and obtain the output Φ. In particular, the power field Φ is constrained to be dependent on a set of "general" parameters that indicate the stage of the reactor’s life cycle: $\mu :=\left(S_t\mathrm{,}{B}_u\mathrm{,}{P}_w\mathrm{,}{T}_{in}\right)\mathrm{.}$ (4.20) • St: The control rod insertion steps, ranging from 0 to 615, represent the movement of the compensation rods from all rod clusters being out (ARO) to fully inserted. This range takes into account the overlap steps.

• Bu: The average burnup of the fuel in the entire core indicates the amount of energy extracted from the fuel and increases over time. Its value ranges between 0 (at the beginning of the fuel’s life cycle) and B_{u, max}, which is set to 2500 MWd/tU (the end of the fuel’s life cycle). The specific evolution of Bu depends on the reactor’s operational history.

• Pw: The power level of the reactor core ranges from 0.3 to 1 FP (full power).

• Tin: The temperature of the core coolant at the inlet falls within the range of 290 to 300°C.

In this manner, Φ can be implicitly represented by $\varphi \left(\mu \right)=\varphi \left(S_t\mathrm{,}{B}_u\mathrm{,}{P}_w\mathrm{,}{T}_{in}\right)$, thanks to CORCA-3D. There are 177 fuel assemblies in HPR1000, and each assembly is numerically represented using 28 vertical levels. Thus, Φ is a vector with dimensions N=4956(=177 × 28). The discrete solution set 𝑴=$\left\{\varphi \left(\mu \right)\in {ℝ}^N|\mu \in D\right\}$ consists of P=18480 solution snapshots with the parameter configuration $D:=B_u^s\otimes S_t^s\otimes P_w^s\otimes T_{in}^s$, where $S_t^s=\left\{\mathrm{0,}\mathrm{1,}\dots \mathrm{,}615\right\}$, $B_u^s=\left\{\mathrm{0,}\mathrm{50,}\mathrm{100,}\mathrm{150,}\mathrm{200,}\mathrm{500,}\mathrm{1000,}\mathrm{1500,}\mathrm{2000,}2500\right\}$, $P_w^s=RU3\left(\mathrm{30,}100\right)$ and $T_{in}^s=RU3\left(\mathrm{290,}300\right)$. The operator RU3(a, b) represents three independent and identically uniformly distributed samplings in the closed set [a, b] where a<b. In this study, we selected 90% of 18480 snapshots for training the forward and inverse models, and the remaining 10% of snapshots were used for test purposes.

Setup of observation data and modeling the observation noise. In the context of RODT, the observations denoted by Y_o are used to infer the input parameter vector μ and subsequently determine the resulting power field Φ. The arrangement of these observations, organized node by node, is shown in Fig. 6. Each observation value represents the average value over the corresponding node in which the SPNDs are located. It is important to clarify that the observations used in the analysis were not obtained from direct engineering measurements but were derived from numerical simulations of Φ using CORCA-3D.

The noise associated with the observations can be modeled as either a Gaussian or a uniform distribution, depending on the specific requirements. Specifically, the noise term can be mathematically represented as follows: $e=Y_{\text{o}}\cdot N\left(\mathrm{0,}\sigma \mathrm{,}{Y}_{\text{o}}\right)\text{or }{Y}_{\text{o}}\cdot U\left(-\delta \mathrm{,}\delta \mathrm{,}{Y}_{\text{o}}\right)$ In this context, the vector N(0, σ, Y_o) denotes a collection of random variables following a normal distribution. These random variables have the same dimensions as Yo and are characterized by a mean of zero and standard deviation of σ. Similarly, the vector $U\left(-\delta \mathrm{,}\delta \mathrm{,}{Y}_{\text{o}}\right)$ represents a set of random variables uniformly distributed within the interval [-δ, δ]. These variables have the same dimensions as Y_o.

Considering the prevalence of Gaussian noise in real-world nuclear reactors, our subsequent analysis primarily focused on evaluating the performance of the algorithm under Gaussian noise conditions.

Software implementation We integrated our software into the RODT platform using the following modules:

• rodtPro: This module implements data cleaning and preprocessing techniques based on domain expertise.

• rodtROM: This module constructs ROM through techniques like POD, SVD, and RB.

• rodtML: This module focuses on surrogate modeling and inverse model optimization. It also includes hyperparameter optimization techniques.

• rodtPost: This is a post-evaluation module that defines convergence metrics and conducts comparative analyses using evaluation methodologies derived from peer-reviewed literature.

• rodtUI: This module is responsible for visualizing the results, providing an intrinsic review of the methods, showcasing the optimized hyperparameters, and designing user-specific visualizations.

Practice Metrics for Inverse Model Evaluation

To evaluate the performance of the algorithms, we designed a series of numerical experiments using nuclear reactor data produced by an intrusive program with added noise. We carefully selected metrics related to accuracy and finally defined the normalized reconstruction fields L₂ error (E₂) and L_∞ error (E_∞) to evaluate the performance of the comparative algorithms.

The reconstructed field L₂ error is defined in Eq. 4.21 to evaluate the average estimation of the field, where ${\Vert •\Vert }_2:={\left({\displaystyle {\sum }_i}{\left({•}_i\right)}^2\right)}^{1/2}$ denotes the L₂norm. $E_2:={\Vert Q{F}^{\text{ML}}\left({\mu }^{*}\right)-{\Phi }_{\text{true}}\Vert }_2/{\Vert {\Phi }_{\text{true}}\Vert }_2$ (4.21) The reconstructed field L_∞ error is defined in Equation 22 to evaluate the worst field estimation, where $\parallel •{\parallel }_{\infty }:={{\displaystyle \max }}_i\left|{•}_i\right|$ denotes the L_∞norm. $E_{\infty }:={\Vert Q{F}^{\text{ML}}\left({\mu }^{*}\right)-{\Phi }_{\text{true}}\Vert }_{\infty }/{\Vert {\Phi }_{\text{true}}\Vert }_{\infty }$ (4.22) It should also be noted that the reconstructed observation L₂ error and L_∞ can be generated similarly.

Accuracy and Stability Analysis of the Parameter Identification Phase

In this section, we describe the tests conducted to evaluate the accuracy of parameter identification. Initially, we assessed the effectiveness of the proposed hybrid methodology, the coarse–fine–grid search, as discussed in Sect. 3.6, using the KNNLHS method as an illustrative example. To visualize this process, we generated scatter plots depicting the predicted values of the parameters St and Pw, which were identified as the two dominant parameters in a previous study [15] (see Fig. 7 for six cases with different input parameter settings.

Fig. 7Full-screen View

Distribution of generated candidate parameters ${\mu }^1=S_t$ and ${\mu }^3=P_w$ with hybrid algorithms. The solid green points are the true values, while the hollow yellow points stand for the candidate parameter generated by the mere KNN algorithm. The hollow blue triangles represent the five optimal samplings by LHS around the initial guess given by the KNN evaluated in the observation space. The solid red points represent the output of the algorithms, which is the best among the five optimal samplings

Figure 7 shows a visualization of the local search space in the inverse algorithms. The figures in the first row depict scenarios with varying Pw values ranging from relatively low to moderate to relatively high. Correspondingly, the figures in the second row illustrate cases with varying St values, again encompassing relatively low, moderate, and high values. These figures effectively show the considerable improvements achieved using the KNNLHS method, particularly in the simultaneous optimization of both Pw and St. The visualizations underscore the effectiveness of the method in yielding enhanced optimization outcomes.

The parameter prediction deviation is demonstrated in Fig. 8 over the test set. Owing to space limitations, we chose only KNN, KNNLHS, and FSA for illustration and tested clean and noisy (σ= 5%) observation data. Using the LHS algorithm, the KNN prediction for parameter identification became more accurate, and the prediction deviation for the FSA algorithm was concentrated at relatively low intervals, thus proving its excellent accuracy. In general, the accuracy of burnup, the power level, and the control-rod step are acceptable from an engineering perspective.

Fig. 8Full-screen View

Prediction deviation of μ with inverse model based on the clean or noised (σ= 5%) observation data

Comparison of the Convergence Rate among Metaheuristic Algorithms

The convergence rate is a crucial criterion for evaluating algorithms that are integrated with the generation design; thus, in this section, we compare different metaheuristic algorithms and their hybrids. Their convergence performance in 50 iterations (generation) with clean and noisy (σ=1% and σ=5%) observation data is shown in Fig. 9.

Fig. 9Full-screen View

Convergence rate of the cost function (defined by the L₂ norm) with different metaheuristic optimizers. The x-axis represents the iteration (generation), and we set the maximum iteration (generation) to 50

When the observation data were clean, the FSA exhibited the fastest convergence rate, followed by the advanced differential evolution (ADE). Moreover, when integrated within an SVC coarse-grid search, both the DE and CSDE algorithms demonstrated the ability to escape the local optima in the early stages while incurring a relatively low-cost function, thereby surpassing their standalone performance.

Among the two population-based algorithms, DE outperformed CSDE. This can be attributed to two factors: (i) the adoption of the best mutation method for DE, where we use the optimal mutation strategy for DE while using the classical strategy; and (ii) the relatively straightforward nature of our problem, wherein the inverse problem lacks multiple optima. Consequently, the CS step in each generation of CSDE may lead to unnecessary exploration of the parameter space that is distant from the optimal solution.

In the presence of noise-contaminated observation data, FSA continues to exhibit superior convergence capabilities. However, its advantage diminishes as the level of Gaussian noise reaches 5%, eventually leading to a decline in advantage. Notably, the ADE and advanced cuckoo search differential evolution (ACSDE) failed to deliver satisfactory results under noisy conditions. This is primarily due to the erroneous predictions of SVC, which further leads to inadequate initial estimates provided by population-based algorithms.

Evaluation of State Estimation Phase by Accuracy

In this section, we evaluate various methods in the state-estimation phase based on accuracy. First, we compared the accuracy and stability under two metrics: the relative error in the L₂ norm defined by Eq. 21 and the L∞ norm defined by Eq. 22 for field Φ, which is presented in Table 2. The ‘Mean’ values in the tables represent the average accuracy across 40 tested specimens, while the ‘STD’ values indicate the standard deviation of errors over the same set of specimens. The optimal values are indicated in bold.

In general, integration of the LHS algorithm enhanced the accuracy of the KNN model. This improvement can be attributed to the LHS sampling approach, which mitigates the discontinuity issue associated with KNN predictions by sampling around the initial estimate.

Furthermore, the NN method exhibited a decline in accuracy when confronted with noise-contaminated observational backgrounds. This is attributed to its reliance on gradients, rendering it less effective in discrete cases.

To visualize the performance of the inverse algorithms directly, we randomly selected one specimen from a test dataset. Figures (11 and 12) depict the relative L₂ error of the field prediction at the 11th vertical level. Readers are referred to Fig. 10 for the true field values.

Fig. 11Full-screen View

(Color online) Reconstructed relative errors of radial power distribution on the 11th axial plane with clean observations

Fig. 12Full-screen View

(Color online) Reconstructed relative errors of radial power distribution on the 11th axial plane with observations of noise level σ= 5%

Fig. 10Full-screen View

(Color online) True value of the power distribution on one specimen in the test dataset over the core in a realistic HPR1000 reactor

To fully demonstrate the reconstruction accuracy of our approach, we examined cases in which the reconstruction was relatively poor, particularly around the control rod interfaces. Reconstruction of the physical fields in locations near the control rod interfaces is challenging because the behavior in these regions is not easily captured with high fidelity. Readers can refer to Fig. 10 for the true field value on this axis, and Figs. 13 and 14 depict the prediction errors of various methods on this axis.

Fig. 13Full-screen View

Predicted relative error for the D01 fuel assembly with clean observations

Fig. 14Full-screen View

Predicted relative error for the D01 fuel assembly with observations of noise level σ= 5%

However, our results also show that even in such difficult cases, our method can achieve good reconstruction in most cases. Although some techniques performed worse around the interfaces, the majority of our proposed methods still facilitated a reconstruction quality comparable to other areas across the domain. Overall, this validation indicates that the proposed approaches hold promise for representing complex reactor physics behaviors with suitable accuracy, even when challenges arise locally because of design complexities, such as control rod insertion points.

Evaluation of State Estimation Phase by Time Cost

In this section, we present an overview of the time costs (we take the mean of the 40 specimens) associated with each proposed algorithm in the state-estimation phase, which can be found in Table 3. The optimal values are highlighted in bold, and the second optimal value is underlined.

It is worth mentioning that the NN method demonstrated notable speed. However, metaheuristic algorithms incur relatively high time costs. This is primarily attributed to the repeated computation of the cost function in Eq. 12, which occurs 50 times for each population member over 50 iterations/generation cycles, resulting in a total of 2500 cost function evaluations. Moreover, these algorithms involve additional computations and sorting operations for information exchange and identification of the best solutions, contributing to the discrepancy in time cost.

Furthermore, it is important to note that the methods LHS, ALHS, and LHS2STEPS have approximately the same time costs. In the algorithm LHS, we set the number of samples to 1000. In algorithms ALHS and LHS2STEPS, we set the first-stage sampling number to 750, and the second-stage sampling is repeated five times, with each repetition sampling number uniformly set to 50, adding up to 1000 in total. Model-based algorithms, such as KNN and its hybrid KNNLHS, have relatively low time costs, that is, less than one second, which is acceptable in the industry.

Comprehensive Evaluation of Various Algorithms

In the state-estimation phase, a comparison of the proposed algorithms in terms of time cost and accuracy is presented in Fig. 15. The KNNLHS algorithm demonstrates a commendable level of accuracy, achieving solutions within 0.1 s for both clean and noise-contaminated observation backgrounds. Consequently, KNNLHS can be considered suitable for real-time online inverse modeling in the context of RODT. Moreover, in scenarios in which time constraints are not limiting factors, the FSA, DE, and ADE algorithms provide reasonably accurate solutions.

Fig. 15Full-screen View

Comparison of the synthetic performance for different inverse algorithms