Automated resonance evaluation; non-convex decomposition method for resonance regression and uncertainty quantification

. This work serves as a proof of concept for an automated tool to assist in the evaluation of experimental neutron cross section data in the resolved resonance range. The resonance characterization problem is posed as a mixed integer nonlinear program (MINLP). Since the number of resonances present is unknown, the model must be able to be determine the number of parameters to properly characterize the cross section curve as well as calculate the appropriate values for those parameters. Due to the size of the problem and the nonconvex nature of the parameterization, the optimization formulation is too difficult to solve as whole. A novel method is developed to decompose the problem into smaller, solvable windows and then stitch them back together via parameter-cardinality and parameter-value agreement routines in order to achieve a global solution. A version of quantile regression is used to provide an uncertainty estimate on the suggested cross section that is appropriate with respect to the experimental data. The results demonstrate the model's ability to find the proper number of resonances, appropriate average values for the parameters, and an uncertainty estimation that is directly reflective of the experimental conditions. The use of synthetic data allows access to the solution, this is leveraged to build-up performance statistics and map the uncertainty driven by the experimental data to an uncertainty on the true cross section.


Introduction
Resolved resonance range (RRR) cross section evaluation is a process that takes a long time partly due to resonance identification being performed manually by expert nuclear data evaluators. The current practice of a generalized linear least squares (GLLS) fitting procedure does not consider discrete variables, making this part of the fitting procedure laborious, timeconsuming, and irreproducible from evaluator to evaluator. Additionally, it is well-documented that the uncertainty estimation given by the default GLLS procedure is unreliably low [1]. Having these uncertainties that are "too small" results in additional work for the resonance evaluator to search for a nonsystematic, artisanal, way to inflate the estimated uncertainty to be more appropriate. The goal of this work is to step through a proof-of-concept for creating a tool that addresses these issues. Such a tool must be able to automatically determine 1) the number of resonances present, 2) the discrete spin group assignment for each, 3) mean values for the continuous resonance parameters that describe each of the identified resonances, and 4) a measure of uncertainty on the resonance parameters that is reflective of the experimental data uncertainty.
The capture cross section of the copper-63 isotope is the focus for this demonstration. The methodology is developed at a decreased level of physical fidelity, namely, the Single Level Breit-Wigner approximation * Corresponding author: nwalton1@vols.utk.edu to the cross-section parameterization is used and only a single spin group is considered. As discussed later in this article, the ability to synthesize data allows for both the tool and the datasets to reflect these physical approximations. The following methodology and results demonstrate the potential to integrate modern computation methods/machine learning (ML) methods into the resolved resonance range evaluation process. The motivation is for these methods to serve as a fast, accurate, and reproducible tool for RRR evaluators to have a starting place in their work. The work presented here is based on the master's thesis in reference [2].

Data synthesis
For exploratory analysis and development of any methodology, it is useful to have a large, labelled data set. Neither of which is the case for experimental neutron cross section data. There are a limited number of RRR experimental data sets, none of which can be considered "labelled". Labelled data indicates that the solution is accessible, were that the case for neutroncross sections, any further measurement or evaluation would be unnecessary. This challenge can be addressed by generating synthetic datasets, a common practice in many ML and statistical applications.
For this problem, the physical model is the R-Matrix theory for neutron resonance reactions [3]. This is a well-established theory, and the input parameters are observed to fall on well characterized statistical distributions. The approach to data synthesis then becomes: 1) sample a realization of the resonance parameters and take them as determined or "true", 2) calculate a cross section from the true parameters, 3) add experimental noise. This process produces a noisy experimental dataset with a corresponding solution.
Leveraging this data synthesis method allows for control over the data complexity. In this case, data was generated using approximations to the R-Matrix theory, allowing for a much more achievable proof-of-concept challenge. The exact fidelity of the synthetic data is as follows. The Porter-Thomas and Wigner distributions were used to sample resonance parameters [4]. These distributions describe values for the resonance widths and spacings relative to the average resonance width/spacing of the nucleus. These average values for the copper-63 isotope were taken from Mughabghab [4]. Only the first s-wave spin group was considered with a total orbital angular momentum of 1 and negative parity. To calculate the theoretical cross section, the sampled parameters, spin group, and isotope specific constants were input to the Single Level Breit-Wigner (SLBW) approximation to R-Matrix theory [1]. This approximation greatly simplifies the functional form but maintains the basic structure of resonances. In practice, several experimental corrections to the theoretical cross section must be made to produce a cross section that could be seen in a laboratory setting (i.e., background signal, doppler broadening, etc.). For this work, these corrections are neglected.
Due to several experimental effects and the nature of the measurements themselves, a varying amount of noise/uncertainty is seen across the energy dependent cross section. To approximate this, a normal distribution with a dependent mean/standard deviation was used to sample a noisy experimental data vector around the theoretical cross section vector: Where !"#$ is the theoretical vector and #%& is the experimental vector. Similarly, the uncertainty assigned to the experimental data points is given by the dependent standard deviation term. The a and b parameters control the level of noise/uncertainty. This distribution was selected as an approximation to the noise model for counting statistics, the Poisson distribution, which resembles ( , √ ) at or above approximately 30 counts. While in practice, the cross section is not measured as counts directly, the value-dependent noise/uncertainty is representative of the capture measurement where the areas of higher cross section values (resonance peaks) will see smaller relative noise and uncertainty than areas of lower cross section values.

Regression and window decomposition
The resonance identification and characterization problem was framed as a constrained, mixed-integer, non-linear program (MINLP). The objective being to minimize the mean squared error between the calculated cross section and the experimental data points, this is shown below in equation 2.
The subscript i indexes the independent axis in units of energy, the experimental data at energy ' is denoted as ()* ! , and the SLBW calculated cross section is denoted as ) ' . The parameters + , Γ , , -. are inputs to the SLBW cross section formulation and represent the energy location, the capture width, and the square of neutron reduced width amplitude of a resonance respectively. Finding these parameters is the goal of the regression. As mentioned, only a single, s-wave spin group is present in the synthesized experimental data and therefore only that spin group is used in the regression problem. The SLBW formalism expresses the cross section in some energy domain as a simple sum over each of the resonances. Because of this, the input parameters can be thought of as vectors with each index corresponding to a resonance.
' denotes the energy at which the cross section is being evaluated and k represents the index for each resonance present in the domain of interest. The following constraints were applied to this problem.
. , < Γ *,, < . , (,, < (,,53 The L and U values give the lower and upper bounds to each parameter while the / values represent a binary variable that is 1 if the resonance is present or 0 if the resonance is not. This allows for a semi-continuous constraint in the following way: . Where Res Par, L, and U are stand in variables for any of the resonance parameters and their respective lower or upper bounds. The idea is to solve this problem with parameter vectors that are known to be longer than necessary, then the binary variable ( / ) will eliminate a number of those variables depending on the predicted number of resonances within the domain. The L and U bounds are physics informed. The width parameters are bounded by the limits of the Porter Thomas distribution and the energy location parameter is bounded by the energy domain being evaluated. Finally, the constraint described in equation 7 enforces that the index of resonances (k) must be ordered with respect to the energy location of that resonance. Now that we have established a proper formulation of the optimization problem and its constraints, we must input to a solver that can handle it. The solver of choice was Baron [5], a global optimization solver that leverages a powerful branch and reduce algorithm that can accept the nonlinear, mixed-integer constraints as is required by this problem formulation.
The regression problem outlined above describes the general objective function being solved. Solving this can still be challenging even for state-of-the art optimization software such as Baron. Resonance region measurements often contain thousands of data points in energy space and hundreds to thousands of resonances, making a single solve intractable. To address this challenge, a novel method was developed to decompose the energy domain into multiple, overlapping windows that are solved independently and then stitched back together. The overlapping region helps to prevent overfitting in any one window as all windows are forced to agree on a solution to the overlapping regions [2].
For the independent windows, the formulation described in section 2.2 is simply updated to reflect energy constraints that correspond to the determined window domain. Once each window has arrived at a solution, the overlapping regions are forced to agree on a single, global solution. This is done in two steps, 1) a cardinality agreement routine that enforces agreement on the number of parameters used to fit a region (analogous to the number of resonances identified), then 2) a parameter agreement routine that enforces agreement on the values for each parameter. In both cases this is done through the introduction of a penalty term and the constraints remain the same. The exact formulation can be found in reference [2].

Uncertainty estimation
The optimum solution to the problem formulation in section 2.2 is a predicted number of resonances and respective parameters. Because of the nature of experimental measurements, the true underlying cross section is unknown and therefore any evaluation must also give an estimate of uncertainty on the suggested parameters. These uncertainties are vital to the community's confidence in simulations of nuclear systems and can have a significant impact on the safety, efficiency, and engineerability of any nuclear system. The methodology developed in this article estimates uncertainty through a quantile regression routine. This is a well-established method that expands upon linear least squares regression by estimating the median, or other quantiles, rather than the mean and does not require conditions such as linearity, normality, or constant variance [6]. The quantile regression minimizes the median absolute deviation across a dataset for a given quantile rather than the mean squared error. The 0 function is defined as a "check" function that gives an appropriate weight for the quantile of interest ( ) based on the error in the data points (u). The inverse standard deviation of the experimental data points is introduced to help with nonlinearities [2]. The different quantiles correspond to the number of experimental data points that are over predicted, i.e., the 10 th quantile would overpredict 10% of the experimental data points and underpredict 90%, likewise, the 50 th quantile corresponds to the median.
Solving for different quantiles can inform a confidence envelope around the mean value and variances on the resonance parameters, informed by the experimental data points and associated uncertainties. This method differs from the current practice of GLLS in its ability to handle non-linearity, asymmetric distributions, and changes in the variance.

Data synthesis and regression
Large datasets with many resonances were generated using the data synthesis methodology in section 2.1, then the space was separated into independently solved windows. The selected optimization solver was able to handle the MINLP well, often predicting a cross section with a smaller mean squared error than that of the solution. An example of this is shown in figure 1. The regression routine also saw significant success identifying resonances in seemingly unstructured data, as shown in figure 2. Because the solver starts with a given number of resonances and then eliminates them if necessary, the maximum number must be decided on. Through development and implementation of this methodology, the best performance was seen with a maximum of 5 resonances per window. Extreme limits of the statistical distribution for resonance spacing were used to define energy windows with confidence that no more than 5 resonances would be present.

Uncertainty estimation
The primary goal of the quantile regression routine is to estimate an uncertainty on the predicted cross section that is representative of the noise and uncertainty in the experimental data by solving for quantiles that over or under predict different percentages of the experimental data points. Upon further inspection of the quantile solutions with respect to the experimental data, a mapping was found to the quantiles with respect to the true cross section. Figure 3 plots the ratio of overpredicted true cross section values (solution quantiles) with respect to the experimental quantiles at each value of . The visually identifiable functional form is that of the logistic function. This result means that an uncertainty envelope around the experimental data can be mapped to an uncertainty envelop around the true cross section and is done so with a well-characterized statistical distribution. This claim, of course, depends on the fidelity of the synthetic data.
Beyond the quantile regression routine, unlimited access to data/solution pairs presents the opportunity for several other supervised machine learning or statistical methods. One simple example is identifying weaknesses or bias in the regression routine through the build-up of performance statistics. This is shown in Figure 4.

Window stitching routine
The window decomposition method proved to be a reliable way to make this a tractable optimization problem. The recombining of windows, or stitching, through cardinality and parameter-value agreement routines saw success in several cases. For example, the iteration was able to arrive at proper solutions for cases where the two overlapping windows disagreed on the number of resonances present, were both incorrect about parameter cardinality or parameter values, or even when resonance structures were cut-off by the energy boundaries of a given window [2]. The data synthesis technique allowed for targeted testing of each of these specific cases. This iterative routine was found to be reliable because its success was driven mostly by the ability to accurately solve the independent windows. If the window solutions are tractable, this routine accurately produces an updated, global solution.

Conclusion
This work highlights the results of a promising demonstration for the development of an automated tool to identify and characterize resolved resonances in experimental cross section data. This would reduce the manual workload required and introduce reproducibility into the evaluation process. The MINLP formulation allows for the number of resonances present to be considered, while the method for window decomposition makes the optimization problem tractable. The method of data synthesis provides two major advantages, 1) control over the domain for exploratory analysis/development, and 2) the opportunity for statistical or ML methods that have the potential to accurately characterize the uncertainty on the true underlying cross section. While this methodology was carried out with lower-fidelity physics than would be used in practice, the results motivate an extension of this work to a higher-fidelity evaluation.
This material is based upon work supported by the Department of Energy National Nuclear Security Administration through the Nuclear Science and Security Consortium under Award Number(s) DE-NA0003180 an d/or DE-NA0000979.
This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States