MUF and SigmaMUF#
Historical context#
Fulfilling safeguards regulations and agreements requires demonstrating that nuclear material has not been lost, removed, or otherwise been unaccounted. Both containment and surveillance methods in addition to direct accountancy of nuclear material is used to meet these requirements. Containment and surveillance (C/S), as the name implies, is used to contain (e.g., seal containers) and surveil (e.g., optical cameras) nuclear material. C/S is complemented by quantitative material accountancy which seeks to quantify the amount and form of nuclear material in a given area. MAPIT focuses on providing tools and routines used in accountancy of nuclear material. The material balance, discussed in depth here, is the cornerstone of nuclear material accountancy. For a longer overview of the goals and safeguards and C/S or the history of material accountancy, readers are encouraged to check the references as a detailed history of these topics are out of scope for this document.
Theory#
MUF#
Accountancy of nuclear material is of interest to many regulatory bodies. One principle quantity used to ensure material is accounted for is the material balance calculation. This is sometimes also called Material Unaccounted For (MUF) or Inventory Difference (ID). Material balance calculations are performed over defined physical areas of a nuclear facility, called material balance areas (MBAs), at regular intervals called material balance periods (MBPs). There are a variety of metrics and criteria used to determine both the MBA (in both structure and quantity) for a given facility and the associated material balance period. Discussion of MBA selection criteria will not be discussed here and we refer the reader to the references for further inquiry. The material balance can be calculated by understanding the inputs, inventories, and outputs for a given material balance area. First, let the material balance period be represented as a non-zero, positive real integer:
Next, consider the sequence of material balance period values:
Then the
Alternatively, if continuous flows are present (i.e., continuous inputs and outputs), then the material balance can be represented as:
is the input to the material balance area at time and location is the output of the material balance area at time and location is the inventory at time and location was chosen to denote container and avoid overloaded notation between inventory and input terms
The material balance for has three primary terms; input, output, and inventory. For each term, there may be multiple different locations so
There is a separate expression for material balances with discrete items versus continuous flows. The expression for the latter uses integrals over the time period of interest to denote that these quantities are usually integral rather than summed. For example, discrete canisters of material arriving to a material balance area might have their contents weighed, assayed, and summed for the material balance period. A continuous input might be a flow measured in mass per unit time, which would be integrated over the material balance period instead of summed.
Note
Material balances are calculated at discrete intervals of time and are consequently not continuous. There are a variety of ways to represent the material balance graphically. In MAPIT, we opt to have a continuous representation by holding the value calculated at


The material balance is intuitive and does not make assumptions about potential material loss pathways. The material balance should be exactly zero due under normal operating conditions as all material should be accounted for. In contrast, the material balance should be non-zero under anomalous conditions that cause material losses or gains that are not measured. However, the material balance will always be non-zero, even under normal operating conditions, due to the presence of measurement error. Additional analyses are then required to detect material loss in the presence of measurement uncertainty.
It is useful to express a series of MUF values as a sequence which facilitates various trend testing and statistical analyses:
The MUF sequence can be represented as a Gaussian distribution as, over long enough periods of time and enough samples, even measurement biases behave as random errors. The fundamental goal of material accountancy then is to detect a shift in the distribution of values in a MUF sequence. Consider the two normal distributions below with different means and a standard deviation of one. This shift in the mean between two distributions represents behavior that would occur during an anomalous operation of a facility. However, the mean shift between the distribution would be difficult given that the shift is small compared to the standard deviation of the distributions.
Tip
You can verify this yourself by calculating an odds ratio using Bayesian principles. That is, for two normal distributions
See also
The guided exercises included in the documentation contain concrete examples of the impact of measurement error on the MUF distribution and overall detection probability.


Now consider the same mean shift, but with a smaller standard deviation for both distributions. The overlap between the distributions has become notable smaller which makes discriminating between the two distributions easier. Any loss pattern, and thus the underlying anomalous MUF distribution, would be difficult to quantify in practice. However, this example should make clear that the uncertainty of the MUF sequence is an important factor in detecting material loss. There are a variety of techniques that can be applied to the MUF sequence to monitor for anomalous behavior, but all techniques are ultimately limited by the uncertainty.


MUF#
Note
For simplicity, the following derivation shows the calculation of a single entry of the
For simplicity we introduce the notation of
Substituting in the multiplicative error model for each term and starting with the input:
Note the following:
behaves as a constant for a specific slice in timeThe true input won’t change at an instant in time regardless of how many times it is measured
Continuing on and zeroing out constant terms it follows that
It’s generally assumed that the random and systematic errors are random variables, and consequently uncorrelated with each other, which leads to
As the true input,
Assuming that each input location is uncorrelated, then it is possible to simply sum the variances for each location together such that
A similar exercise can be performed on the output term leading to the expression
The inventory term differs from the inputs and outputs since there is a temporal correlation between inventory
It is assumed that the constant values for inventories are uncorrelated, random and systematic errors are uncorrelated with each other, and random errors from different times are uncorrelated. However, systematic errors at the same location but different time, (i.e.,
This is again approximated using the actual measured values as the true values are unobservable and is summed across locations.
Substituting expressions for variance of input (5), inventory (7), and output (6) into the definition of
Note
This assumes one strata for each measurement. That is, an item or flow is measured once when it has the same attributes. Performing multiple measurements on the same strata will reduce the relative standard deviation terms by approximately
Equation (8) can be expanded by replacing the inputs and outputs with their non-integrated quantities. This is done to better illustrate the different components that are calculated in MAPIT.
Discussion#
Much of the traditional literature and research around material balances and associated testing was developed in the 1980s. The relevant seminal papers are listed in the additional resources portion of this document. There are a few key points from these papers that are important to note.
A natural inclination to improve performance of any testing on the material balance would be to reduce the overall uncertainty through 1) smaller material balance areas which result in smaller inventory terms, 2) more frequent material balance frequency resulting in smaller input and output terms, or 3) some combination of the two. Avenhaus and Jaech considered this question and found that none of those procedures necessarily lead to better performance, and in some instances, might lead to a lower detection sensitivity. The work by Avenhaus and Jaech was particularly notable as it lead to several important findings:
Considering a statistical test with maximum test power (i.e., a test with the highest probability of detection for any loss of material of a particular size) applied to a fix length of time, the optimal test is one that ignores intermediate balance evaluations. Put another way, the optimal test for a loss of material over a fixed period of time is one that considers the sum of all intermediate MUF values. This is the same as if no intermediate MUF values had been taken; as if the balance was conducted over the entire period of interest. This applies to protracted, but not abrupt, losses. Here, protracted loss is one that occurs over multiple balance periods or areas.
Concrete example; there might be a regulatory goal for detecting a loss within 3 months. The optimal statistical test would be a balance over a 3 month period; performing a monthly balance provides no benefit with respect to maximum detection probability of a loss over a three month period.
Combining intermediate MUF values in some optimal way, perhaps as a weighted average, still results in a lower detection probability than a global MUF that only considers the beginning and ending states. * An important exception here is that this statement only applies to the unknown loss pattern. Performance improvements can be seen with an optimally ordered MUF sequence if the loss pattern is known, but in practice, the loss pattern is never known
Even applying statistical tests to each intermediate MUF value and then linearly combining the results (as opposed to a test on the combined MUF values) still results in a lower detection probability than a test on a global MUF.
Further, Avenhaus and Jaech showed that not only do these statements apply to time (i.e., different material balance periods), but also to space. Subdividing a material balance area provides no benefit in terms of detection probability, and perhaps even decrease detection probability, compared a test on the larger material balance. These statements also assume a fixed false alarm probability.
Important
Avenhaus and Jaech’s work applied only to the probability of detection. There are other benefits to subdividing material balances into smaller units of time or space; principally to localize a potential material loss in space or time, but this comes at a cost to detection probability, specifically for protracted losses that are split over multiple balances or areas. Again, this statement applies only to protracted losses, not abrupt. Avenhaus and Jaech’s work only considered random errors, but Burr and Hamada later went on to show that the inclusion of systematic error does not change the limitations of subdividing material balance areas (i.e., probability of detection does not increase for more frequent balances and smaller balance areas).
Material balance iterations#
In practice, only one material balance sequence,
MUF implementation#
Important
Note that this function is not intended to be used standalone through direct calls, rather, it is designed to be called through the MBArea
class of StatsProcessor
.
The MAPIT implementation of MUF is the first function in MAPIT.core.StatsTests
. The MUF calculation requires a few key variables:
Input, inventory, and output measurements
Input, inventory, and output times
Together the measurements and times should from a timeseries. The time entries should represent the time at which the measurement is taken and should monotonically increase. So a measurement taken at the start is t=0 and a measurement taken one day later should be t=24 (hours) or t=1440 (minutes), etc. There are no unit requirements, but the timeseries should all use the same units of time.
Material balance period
This should have the same units as the input, inventory, and output time
Since each measurement provided to MAPIT could potentially have a different length and/or number of time steps, we first determine the maximum timestep:
117 #out how long the imported data is
118 #as it cannot be assumed (len(data)) = time
119
120 iterations = inputAppliedError[0].shape[0]
121
122 A1 = np.max(np.asarray(list(chain.from_iterable(processedInputTimes)))) #unroll list
MAPIT calculates the entire MUF and
Example
Consider the following example
MBP = 100
Largest time in dataset: 670
Iterations: 25
MAPIT will calculate a
Important
MAPIT represents periodic statistical quantities as continuous, so although there are only 6 balance periods, this is represented as 600 timesteps (once per unit time). Each material balance iteration (i.e., slice [n, 0:600]) will only have 6 unique values that are held constant between material balance updates.
MAPIT calculates the distribution of MUF values (i.e.,
At a high level, MAPIT vectorizes the error model iterations, but uses a for-loop to iterate over the location and time components.
Of the two loops (balance periods and locations), MAPIT uses the time component as the outer for-loop and locations as the inner for-loop.
Note
The outer time loop is indexed from i
to MBPs
rather than starting at zero. This facilitated more understandable indexing, particularly when slicing some of the time indices, but runs counter to the standard python notation that indices start at 0.
The outer for-loop is defined as follows:
150 #-------------------------------------------------------------------------------#
The outer loop calculates the muf
The summation component represents the inner for-loop and is split across three separate loops, each representing the measurement type:
154
174 processedInputTimes[j] <= MBP * i).reshape((-1,))
192 if ispar == False:
The integral component, assumed to be flows in units of mass/time, need to be integrated before they can used in the balance calculation. This is done using a customized trapezoidal integration routine, MAPIT.core.AuxFunctions.trapSum
, which is explained in more depth in the computational considerations part of this guide.
The relevant segment of time corresponding to logicalInterval
variable (note that the interval for the input terms is shown, but it is also calculated for the output terms):
164 if doTQDM:
165 if ispar == False:
166 pbar.update(1)
The relevant times for the current balance period must also be identified for the inventory, but the procedure is easier as only the start and end points of the time series need to be identified, not all values in the interval. This is because the inventories, already assumed to be in the correct mass units, do not need to be integrated and can be subtracted (i.e.,
Finally, all of the components are combined and MAPIT iterates over material balance periods and locations to calculate MUF.
MUF implementation#
Important
The final expression for
Important
Note that this function is not intended to be used standalone through direct calls, rather, it is designed to be called through the MBArea
class of StatsProcessor
.
The implementation of
The entire sequence from Equation
mufdist
is calculatedInput, inventory, and output terms are considered separately
Time components are vectorized whereas locations and balance periods are expressed in explicit floor loops
The structure of the calculation is identical to that of MUF. The key deference is the quantity calculated, here
Following from the muf calculation, integral terms are integrated using MAPIT.core.AuxFunctions.trapSum
with locations and balance periods iterated over using a for loop. It is often desireable to track the contribution to
Here, the random and systematic contributions can be tracked separately. The input and output terms are calculated in a similar way, but only the input term will be discussed for brevity.
490 pbar = tqdm(desc="Sigma MUF", total=int(totalloops), leave=True, bar_format = "{desc:10}: {percentage:06.2f}% |{bar}| [Elapsed: {elapsed} || Remaining: {remaining}]", ncols=None)
491
492
493 InpVar = np.zeros((iterations,int(MBPs * MBP)))
494 InvVar = np.zeros((iterations,int(MBPs * MBP)))
495 OutVar = np.zeros((iterations,int(MBPs * MBP)))
496
497 #-------------------------------------------------------------------------------#
498 #------------------------------ SEID Calculation -----------------------------#
499 #-------------------------------------------------------------------------------#
500
First, the input term is integrated. Then the integral quantity and user supplied relative standard deviations are squared and summed. These are added together and ‘‘tiled’’. Variables VR
and VS
are a vector with length equal to iterations and must be ‘‘tiled’’ for all time steps in the given balance. The contribution components are then stored in a separate array. This process repeats for all locations and balance periods.
The inventory calculation proceeds in a similar manner but differs in that the first and subsequent balance period calculations differ. The previous inventory during the first balance is assumed to be zero, and rather than trying to account for that in the array by prepending a series of zeros, MAPIT simply drops that term. This is implemented by checking if the first balance period is being calculated (i==1
), and if so, not including inventoryAppliedError[j][:, startIdx]
.
518
519
520 if inputTypes[j] == 'continuous':
521 AFTS = AuxFunctions.trapSum(logicalInterval,processedInputTimes[j],inputAppliedError[j])
522 elif inputTypes[j] == 'discrete':
523 AFTS = inputAppliedError[j][:, logicalInterval].sum(axis=1)
524 else:
525 raise Exception("inputTypes[j] is not 'continuous' or 'discrete'")
526
527 VR = AFTS**2 * ErrorMatrix[j, 0]**2
528 VS = AFTS**2 * ErrorMatrix[j, 1]**2
529