ML14072A211
| ML14072A211 | |
| Person / Time | |
|---|---|
| Site: | South Texas  |
| Issue date: | 03/06/2014 |
| From: | Harrison A South Texas |
| To: | Balwant Singal Division of Operating Reactor Licensing |
| References | |
| MF2400, MF2401 | |
| Download: ML14072A211 (110) | |
Text
NRR-PMDAPEm Resource From: Harrison Albon [awharrison@STPEGS.COM]
Sent: Thursday, March 06, 2014 2:18 PM To: Singal, Balwant
Subject:
STP Sensitivity Report and Presentation for Proposed Meeting Attachments: Sensitivity_Analysis_February_10_2014.pdf; Sensitivity Analysis Framework.pdf
- Balwant, Per our discussion yesterday, here is the background information for a meeting on the STP risk-informed GSI-191 sensitivity analysis. We think that dialog with the staff regarding how the sensitivity was done and the results would be beneficial before we make a formal submittal. We are prepared to meet with the staff at their earliest convenience.
Please call Ernie Kee or me if you have any questions.
- Regards, Wayne Harrison 1
Hearing Identifier: NRR_PMDA Email Number: 1165 Mail Envelope Properties (8C918BCF8596FB49BD20A610FA5920CF01F60589)
Subject:
STP Sensitivity Report and Presentation for Proposed Meeting Sent Date: 3/6/2014 2:17:43 PM Received Date: 3/6/2014 2:19:44 PM From: Harrison Albon Created By: awharrison@STPEGS.COM Recipients:
"Singal, Balwant" <Balwant.Singal@nrc.gov>
Tracking Status: None Post Office: CEXMBX02.CORP.STPEGS.NET Files Size Date & Time MESSAGE 459 3/6/2014 2:19:44 PM Sensitivity_Analysis_February_10_2014.pdf 2146859 Sensitivity Analysis Framework.pdf 511958 Options Priority: Standard Return Notification: No Reply Requested: No Sensitivity: Normal Expiration Date:
Recipients Received:
A Practical Guide to Sensitivity Analysis of a Large-scale Computer Simulation Model
!"
#!"
$ % &"'(")
- +' , &
- & ' , &
. ' , &
-!/ '
. + 0 0', & ,% 1 2 0
& ! ', & ,% 1 2 0
- 3 ' 41 0
A Practical Guide to Sensitivity Analysis of a Large-scale Computer Simulation Model David Morton, Jeremy Tejada, and Alexander Zolan The University of Texas at Austin Abstract We describe a 10-step sensitivity analysis procedure that applies to a large-scale computer sim-ulation model. We propose using tornado diagrams as the initial tool for identifying the input parameters to which the simulations outputs are most sensitive. Sensitivity plots and spider plots complement tornado diagrams by capturing nonlinear responses in outputs to changes in inputs. Regression metamodels, and associated experimental design, help understand sensitivi-ties to, and interactions between, input parameters. Our motivating model from GSI-191 has a number of distinguishing features: (i) The model is large in scale in that it has a high-dimensional vector of inputs; (ii) Some of the models inputs are governed by probability distributions; (iii)
A key output of the model is the probability of system failurea rare event; (iv) The models outputs require estimation by Monte Carlo sampling, including the use of variance reduction techniques associated with rare-event simulation; (v) It is computationally expensive to ob-tain precise estimates of the failure probability; (vi) We seek to propagate key uncertainties on model inputs to obtain distributional characteristics of the models outputs; and, (vii) The overall model involves a loose coupling between a physics-based stochastic simulation sub-model and a logic-based PRA sub-model via multiple initiating events in the latter sub-model. We review a subset of a much larger literature, guided by the need to have a practical approach to sensitivity analysis for a computer simulation model with these characteristics. We illustrate our proposed 10-step procedure on a simple example of a simulation model for system reliabil-ity. Important themes repeat throughout our recommendations, including the use of common random numbers to reduce variability and smooth output analysis, a focus on assessing dier-ences between two model con"gurations, and proper characterization of both sampling error and uncertainties on input parameters. In Appendix A we assess the sensitivity of core damage frequency (CDF) estimates to changes in input parameters for the South Texas Project Electric Generating Station GSI-191 risk-informed resolution project. In particular, we use output from the CASA Grande simulation model to construct a tornado diagram to assess which parameters, from a list of candidate parameters, CDF appears most sensitive, and we further construct a one-way sensitivity plot for one of the most sensitive parameters.
1 Background, Purpose, and Lexicon Paraphrasing Kleijnen [9],
Sensitivity analysis, of a computer simulation model, estimates changes in the models outputs with respect to changes in the models inputs.
In this report, we review approaches to sensitivity analysis of a computer simulation model, fo-cusing on speci"c approaches that we see as practically viable in the context of resolving GSI-191 through a risk-informed approach. And, we propose a 10-step sensitivity analysis procedure. Before proceeding, some remarks on our characterization of sensitivity analysis, adopted from Kleijnen, are in order.
1
- 1. We deliberately use the verb estimates because the true value of our models outputs (e.g.,
probability of system failure) cannot be computed. Rather this output must be estimated using Monte Carlo sampling in our implementation via a computer simulation model.
- 2. For the moment, we are purposefully vague about how to make changes in the models inputs because this is a key distinguisher of dierent approaches to sensitivity analysis. We explore this issue in some detail in this report.
- 3. Another key distinguisher of approaches to sensitivity analysis is the manner in which we quantitatively and qualitatively summarize changes in our estimates of the models outputs.
We discuss some important and complementary approaches.
- 4. We use the plural models inputs because of our interest in being able to handle a high-dimensional vector of input parameters.
- 5. We again use models outputs recognizing that multiple performance measures are of simulta-neous interest. For example, we may be interested in both the core damage frequency (CDF) and the change in core damage frequency relative to a base model (CDF), where the latter model does not account for failure modes associated with GSI-191. Or, we may be interested in the four-tuple (CDF, CDF, LERF, LERF), where LERF denotes large early release frequency.
- 6. There is no redundancy in the term computer simulation model. Physical simulation models dier from models implemented on a computer. The notions of to model and to simulate both involve mimicking a real-world system. However, in our context modeling means the act of abstracting the real-world system into a set of mathematical equations and/or logic constituting an abstract model. While a simulation model is a mathematical model, it is usually taken as distinct from other classes of mathematical models that yield an analytical solution.
- 7. While a simulation model may be deterministic or stochastic, we focus on stochastic simula-tion models. The output of a stochastic simulation model may be deterministic or random.
Consider the probability of system failure in the context of GSI-191. If we condition the input on a random initiating frequency, the output of the model is a conditional probability of system failure; i.e., the output is a random variable. On the other hand, if we integrate with respect to the distribution governing the random initiating frequency, the probability of failure is a deterministic output parameter. We consider both alternatives.
- 8. We cannot compute even a deterministic output measure exactly. Rather, we must estimate the output using Monte Carlo sampling. We have sampling-based errors associated with Monte Carlo methods, but these errors can be quanti"ed. The errors can also be reduced by increasing the sample size, and they can be reduced by using so-called variance reduction 2
techniques. The latter methods are particularly important because of our interest in rare-event simulation in which failure probabilities are very small.
- 9. It is important to distinguish three sources of error: First, we have sampling-based errors associated with Monte Carlo methods discussed above. Second, we have errors due to uncer-tainties on some model inputs. Third, we have errors due to a lack of "delity of the model itself. We can attempt to reduce the second type of error by gathering more data or eliciting (more) information from experts. We can also use data and elicitation to reduce the third type of error if we have competing hypothesized models, or sub-models. In all three cases, sensitivity analysis can help guide where such eorts should focus.
A second key notion regarding sensitivity analysis in the context of decision problems involves understanding which dierences in inputs make a dierence in the decision rather than simply dierences in model outputs. In this context, Clemen and Reilly [3] characterize sensitivity analysis by saying:
Sensitivity analysis answers the question, What makes a dierence in this decision?
. . . Determining what matters and what does not requires incorporating sensitivity anal-ysis throughout the modeling process.
A similar sentiment is re"ected in Eschenbach [6]:
This sensitivity analysis may be used (1) to make better decisions, (2) to decide which data estimates should be re"ned before making a decision, or (3) to focus managerial attention on the most critical elements during implementation.
Should some pipes that currently have "berglass insulation be retro"tted to instead have re"ec-tive insulation to mitigate risk associated with a possible GSI-191 LOCA event? Such an action would incur signi"cant cost and signi"cant radiation exposure to workers. What dierences in input parameters of a model lead to a yes versus no answer to this question? A surrogate question involves the notion of Regions I, II, and III from Regulatory Guide 1.174 [4]. Suppose nominal values of our models input parameters lead to an assessment that the plant is in Region III. Then, we may ask: What changes in the input parameters would lead us to conclude that the plant would instead be in Region I or Region II?
It is not the purpose of this report to directly answer the above pair of questions. Rather, we provide a framework that, when properly applied, can answer these questions. For concreteness, we apply our sensitivity analysis framework to a simple example of a parallel-series system with four components illustrated in Figure 1. Here, our analogous question will be, Should we perform preventive maintenance on component 3 to decrease the probability the system fails prior to a pre-speci"ed time, t0 ? We describe this illustrative example in more detail in Section 2.
3
These are 3 critical components. If 1 2 1/2 either fails, the system fails.
4 Component that may have preventive maintenance Figure 1: The "gure depicts a series-parallel system in which the "rst two components (1 and 2) are in series with a subsystem that has two components (3 and 4) in parallel. One of the two parallel components must be up for the system to be up, along with both components 1 and 2.
2 Illustrative Example Figure 1 depicts a simple model of system reliability we use to illustrate sensitivity analysis in this report. Our emphasis here is on having a concrete model that is rich enough to serve this purpose as opposed to having a high "delity model of an actual system. The block diagram in the "gure depicts four components with independent random failure times T1 , T2 , T3 , and T4 . If a component fails, it will not be repaired. We seek to understand the failure time of the system, given by T = min{T1 , T2 , max{T3 , T4. While T is a random variable, we use two deterministic output measures P{T > t0 } and E[T ], where the former output is our primary performance measure of system reliability; i.e., the probability the system fails after a pre-speci"ed time, t0 . A secondary output is the expected time until the system fails. We have oriented the measures so that we prefer larger values. The parameters of the four random variables, T1 , . . . , T4 , are inputs to our model of system reliability. We assume the four random variables have exponential distributions, and so we have as model inputs the failure rates of each of the components, 1 , 2 , 3 , and 4 (which have units of failures per unit time), along with the time threshold for which we desire the system to survive, which we denote t0 . Usually, we suppress the dependency of T1 , . . . , T4 on their rates but sometimes we write expressions such as T3 (3 ) to emphasize the rate associated with T3 . The inverse of the failure rate is the mean time until a component fails, and often it is more natural to think in terms of these means: 1 1 1 1 1 , 2 , 3 , and 4 , which have units of time. In what follows we interchangeably speak of failure rates or mean times to failure, depending the context. In our example, we have a decision to make for this system. We can operate the system as depicted in Figure 1 with failure rates 1 , . . . , 4 . We call this the base option. Or, at time 0 we can take component 3 o-line, and perform preventive maintenance on that component. Component 3 would then come back on-line at time t = t. Importantly, the system operates even when component 3 is o-line for preventive maintenance. So the system would operate as three 4
components in series with failure rates 1 , 2 , and 4 during the time interval (0, t) and as the system depicted in Figure 1 with failure rates 1 , 2 , 3 /k, and 4 , on time interval (t, t0 ). Here, performing preventive maintenance on component 3 at time 0 reduces its failure rate from 3 to 3 /k for k 1. For brevity we call this latter option the PM option, although it arguably amounts to an upgrade of component 3 given the memoryless property of the exponential random variable. For the PM option, model inputs include 1 , 2 , 3 , 4 , t0 , t, and k. System reliability under the base option is given by: P {min{T1 , T2 , max{T3 (3 ), T4 }} > t0 } . (1) System reliability under the PM option is: P{min{T1 , T2 , T4 } > t}
- P{min{T1 , T2 , max{T3 (3 /k), T4 }} > t0 t}. (2)
We are using the memoryless property of exponential random variables in equation (2), by writing the product and by writing the reliability of the system over time interval [t, t0 ] as the second term. Also, the rates associated with T1 , T2 , and T4 do not change under the base and PM options and hence we do not make them explicit in equations (1) and (2). That said, we investigate below changes in these rates in the course of our sensitivity analysis. We may treat input parameters, such as 1 , as deterministic but vary the parameter for the purpose of understanding the sensitivity of system reliability to changes in 1 . Or, we may treat 1 as a random variable governed, e.g., by a gamma distribution or by a bounded Johnson distribution. In either case, we may compute, or estimate, the conditional output P{T > t0 l 1 } under the base option, where T = min{T1 , T2 , max{T3 , T4 }}. In the latter case, because P{T > t0 l 1 } is a random variable, we may compute, or estimate, the percentiles (e.g., the 5th, 50th, and 95th percentiles) of P{T > t0 l 1 } knowing the corresponding percentiles of 1 . Alternatively, we may integrate with respect to 1 s distribution and obtain P{T > t0 }. Finally, if 1 is governed, e.g., by a bounded Johnson distribution, we could seek to understand the sensitivity of P{T > t0 } to the parameters of the Johnson distribution. In the context of GSI-191, example input parameters include margins governing various failure modes such as the net positive suction head margin for pumps, the structural margin for pump strainers, air intrusion limits for pumps, in-vessel limits on debris penetration, and solubility limits on boron concentration in the core. Other key inputs include the temperature, pH, and water volume of the pool, parameters governing debris generation and debris transport, parameters gov-erning strainer characteristics, and so on. Some parameters such as pool temperature change over time according to a speci"ed input pro"le. In a stochastic simulation model, random variables play a key role and their inputs can be characterized in one of two key ways, and we take the initiating LOCA frequencies as an example. We model a probability distribution as governing the frequency of breaks of various sizes. We can either take as input: (i) the parameters of that probability distribution, which in the case of STPs GSI-191 analysis are the parameters of the Johnson distributions governing the break sizes 5
for the six NUREG-1829 break-size categories or (ii) we can take as input a percentile (e.g., the median) associated with that distribution. This choice aects how we characterize model output. Similar choices can be made for random variables governing, e.g., the strainer safety margins. The discussion above regarding the treatment of P{T > t0 l 1 } is our analog for the illustrative example. In terms of model outputs for GSI-191, we may be interested in both the core damage frequency (CDF) and the change in core damage frequency relative to a base-case model (CDF). Or, we may be interested in the four-tuple (CDF, CDF, LERF, LERF). A detailed physics-based simulation model, such as CASA Grande, can help characterize the risk associated with a speci"c failure mode. However, proper assessment of overall risk requires propagating such failures through a coupled PRA model, and hence proper assessment of sensitivities to changes in underlying input parameters requires a similar propagation. In the remainder of this report we do not discuss a GSI-191 example in detail, even though we are motivated by risk-informed GSI-191 analyses. Rather, we restrict attention to the example discussed in this section to illustrate ideas. This streamlines the discussion and allows us to provide simple and transparent insights on the relative merits of various approaches to sensitivity analysis. 3 A Practical Step-by-Step Guide to Sensitivity Analysis Step 1: De"ne the Model We let f : Rn Rm denote our idealized model of the system. Here, our notation means that the model takes as input the values of n parameters and gives as output m performance measures. The vector of inputs is denoted x = (x1 , x2 , . . . , xn ) and the vector of outputs is denoted y = f (x1 , x2 , . . . , xn ), where y = (y1 , y2 , . . . , ym ). We call the model idealized because we assume, for the moment, that the outputs are known exactly given the values of the estimates; i.e., for the moment we assume we do not need to perform a Monte Carlo simulation in order to estimate the values of the outputs. Our illustrative example has two models rooted in the base and PM options. The base-option model has the following primitives: Four independent exponential random variables, T1 , . . . , T4 , govern the failure times of four components with respective failure rates 1 , . . . , 4 , and the failure time of the system is given by T = min{T1 , T2 , max{T3 , T4 }}. With these constructs and m = 2 outputs we have the model f de"ned by f (1 , . . . , 4 , t0 ) = ( P{T > t0 }, E[T ] ) , where the equations for P{T > t0 } and E[T ] could be further detailed using four-dimensional integrals. (We do not do so here as it does not further our discussion of sensitivity analysis.) An analogous idealized model can be written under the PM option using equation (2). 6
Step 2: Select Outputs of Interest Our model f has m outputs (y1 , . . . , ym ). In step 2 of the proposed process, we can restrict attention to a subset of these outputs. There are a number of possibilities for our illustrative example. We may have m = 1 with (y1 ) = (P{T t0 }) or (y1 ) = (E[T ]) as the single output of interest. Or, we have two outputs of interest: m = 2 and (y1 , y2 ) = (P{T t0 }, E[T ]). We may have m = 3 outputs: (y1 , y2 , y3 ) = (P{T > t0 }, E[T ], P{T > t0 l 1 }), and this can be extended to include additional outputs such as P{T > t0 l i } for all i = 1, 2, . . . , 4. The notion of attribution is tied to our outputs of interest. Consider the base option in our example. Given that our system failed prior to time t0 , we can assess whether this is due to a failure of component 1, component 2, or due to the failure of the parallel subsystem of components 3-4. Thus we can compute P{T1 = T l T < t0 }, P{T2 = T l T < t0 }, and P{max{T3 , T4 } = T l T < t0 }. When focusing on attribution, we could model m = 3 output parameters: (y1 , y2 , y3 ) = (P{T1 = T l T < t0 }, P{T2 = T l T < t0 }, P{max{T3 , T4 } = T l T < t0 }). Of course, we could further assess whether component 3 or 4 caused the failure rather than taking their paired subsystem via P{T3 = T l T < t0 } and P{T4 = T l T < t0 }. Step 3: Select Inputs of Interest Our model f has n inputs (x1 , x2 , . . . , xn ). In step 2 of the process, we have already restricted attention to a subset of the model outputs. It may seem counterintuitive to choose the outputs before choosing the inputs, but this order is purposeful. Our choice of inputs hinges both on what the analyst sees as important and on the outputs of interest that the analyst selected in step 2. The notion of important here is driven by multiple considerations. In our example, the analyst may believe an input parameter may not change whether the base option versus PM option choice leads to higher system reliability until the parameter changes to some relatively extreme value, and the analyst may seek to understand the magnitude of that extreme. Or, the analyst may believe an output depends crucially on an input parameter, and the analyst seeks to understand the direction and magnitude of change in the output with respect to changes in the input. For our examples base option, if we have selected m = 1 with (y1 ) = (E[T ]) then we may choose as the input vector (x1 , . . . , x4 ) = (1 , . . . , 4 ) and drop the time threshold t0 because this is not relevant when estimating E[T ]. If the analyst believes that components 1 and 2 are identical and components 3 and 4 are identical then it may suce to have the smaller dimensional input vector (x1 , x2 ) = (1 , 3 ), because changes in 1 = 2 apply to both components 1 and 2 and changes in 3 = 4 apply to both components 3 and 4. If P{T > t0 } is one of our outputs of interest then we may seek to understand the sensitivity of the failure probability to choices of t0 , and hence include t0 as an input parameter. However, even if P{T > t0 } is one of our outputs of interest we may not seek to understand its sensitivity with respect to t0 , if changes in t0 are highly unlikely or the value of t0 is "xed by mandate. 7
Step 4: Choose Nominal Values and Ranges for Inputs In the previous step, we have selected n input parameters of interest, (x1 , x2 , . . . , xn ). In step 4, we select nominal values for these parameters and lower and upper bounds for each of these input parameters. We denote the nominal values by (x01 , x02 , . . . , x0n ), the lower bounds by (x01 , x02 , . . . , x0n ), and the upper bounds by (x01 , x02 , . . . , x0n ). The nominal value for an input parameter is typically based on the analysts best point estimate for that input. That said, there are sometimes reasons for selecting an appropriately conservative nominal value. Consider our illustrative example. The threshold time, t0 , may denote the lifetime for which we require the system to survive, but we may not know the value of t0 with certainty. We could select a conservative (i.e., large but reasonable) value of t0 , and if P{T > t0 } is suciently close to 1, we may be satis"ed that our system is of suciently high reliability. Sensitivity analysis explores this notion in a richer manner, seeking not just to understand the failure probability at a single, perhaps conservative, value of t0 , but rather to understand the failure probability over a range of values of t0 . Table 1 gives the input parameters associated with our system reliability example. Lower and upper bounds are speci"ed by what the analyst sees as how low or high these parameters might be, in an absolute sense. More typically, ranges are speci"ed so that the interval contains values that are both reasonable and likely (e.g., we might exclude values that have less than a 10% chance of occurring). All seven parameters in Table 1 have absolute lower bounds of 0 and upper bounds of . However, we have no intention of exploring this entire interval. Even if PM might conceivably degrade a component, we will not explore k < 1. Similarly we will not explore t > t0 because under such a large value for the PM time, it is clearly not worthwhile to pursue PM. It is important to choose ranges for the input parameters that the analyst sees as reasonable and commensurate. This task can be dicult, and we do not mean to minimize that diculty. That said, such choices are continually made during the process of modeling a system, and we see this diculty as implicit in the intimate connection between modeling and sensitivity analysis. Step 5: Estimating Model Outputs under Nominal Values of Input Parameters So far we have referred to the idealized model, f (x). So, with m = 1 model output, P{T > t0 }, we can discuss the value of system reliability under the nominal values of the input parameters, x = x0 , given in Table 1. However, for the large-scale stochastic models in which we have interest, we cannot compute f (x0 ) exactly. Rather, we must estimate f (x0 ) using Monte Carlo sampling. Formally this means that we have another model, which we denote, fN (x0 ), where this model is parameterized by a sample size, N . We can compute fN (x0 ), but because its inputs involve sampling random variablessuch as the failure times of the four components in Figure 1the model output, fN (x0 ), is also a random variable. The random sampling error associated with fN (x0 ) can be quanti"ed, provided the sample size, N , is suciently large, using fN (x0 )s standard deviation via the central limit theorem. 8
Table 1: Nominal values, lower bounds, and upper bounds for the input parameters in our system reliability example: 1 1 1 , . . . , 4 denote mean time until failure for the four components depicted in Figure 1; t0 speci"es the desired lifetime of the system; t is the time required to perform PM on component 3; and, if PM is performed component 3s failure rate drops to 3 /k where k 1. Input Nominal Lower Upper parameter (x) value (x0 ) bound (x) bound (x) 1 1 (months) 200 150 250 1 2 (months) 200 150 250 1 3 (months) 50 25 75 1 4 (months) 50 25 75 t0 (months) 18 12 24 t (months) 1 0.5 3 k (unitless) 2 1 5 Using the nominal values for the input parameters from Table 1, estimates of the system reli-ability, expected time to failure, and failure attribution probabilities are reported in Table 2. The table contains point estimates, and estimates of sampling error in the form of 95% con"dence in-terval halfwidths. For example, the point estimate of system reliability is 0.7520 under the base option and 0.7850 under the PM option, as reported in the P{T > t0 } row of Table 2. This suggests that the PM option leads to higher system reliability, however, we cannot ignore sampling error in coming to this conclusion. We also see from Table 2 that the mean lifetime of the system appears to be longer under the PM option. And, we get a sense of how the attribution probabilities change under the base and PM options, with the probability that the parallel subsystem 3-4 is the cause of system failure dropping under the PM option. The three probabilities in rows 3-5 of Table 2 sum to one because they are conditional on the system failing. Under the PM option the attribution to parallel subsystem 3-4 drops, and hence the likelihood of the failure being attributed to components 1 and 2 necessarily grows. Based on Table 2, we are 95% con"dent that the true value for system reliability under the base option lies in the interval (0.7252, 0.7788), and we are similarly con"dent that the true value for system reliability under the PM option lies in the interval (0.7595, 0.8105). We may be tempted to use the two con"dence intervals (0.7252, 0.7788) and (0.7595, 0.8105) to infer that the dierence is not statistically signi"cant (because the con"dence intervals overlap), but this is not the proper way to analyze this dierence. We describe the approach we recommend shortly. 9
Table 2: Estimates of output performance measures for the base-option model and the PM-option model under the nominal values of the input parameters from Table 1. For example, the point estimate for P{T > t0 } under the base option is 0.7520 and a 95% con"dence interval halfwidth is 0.0268. All estimates in the table are based on a sample size of N = 1000. Output Base-option PM-option Measure Model Model P{T > t0 } 0.7520 +/- 0.0268 0.7850 +/- 0.0255 E[T ] 44.0174 +/- 2.3008 54.9770 +/- 2.9976 P{T1 = T l T < t0 } 0.2150 +/- 0.0255 0.2520 +/- 0.0270 P{T2 = T l T < t0 } 0.2350 +/- 0.0263 0.2940 +/- 0.0283 P{max{T3 , T4 } = T l T < t0 } 0.5500 +/- 0.0309 0.4540 +/- 0.0309 It is often signi"cantly easier to estimate the dierence of an output measure under two system con"gurations than it is to estimate the absolute values of that same output measure. When estimating dierences we can take advantage of the simulation technique called common random numbers in which similar components in the two systems see similar inputs. We illustrate this by estimating P{T > t0 } P{TP M > t0 } P{Tbase > t0 } using both common random numbers and independent random numbers with the same sample size N = 1000, and we present the results in Table 3. The table rightly suggests that we can reduce the variance of estimates of such dierences by using common random numbers. Table 3: Estimates of the dierences in output performance measures between the base-option model and the PM-option model under the nominal values of the input parameters from Ta-ble 1 using both common and independent random numbers. For example, the point estimate for P {T > t0 } under the base option is 0.0330 and a 95% con"dence interval halfwidth is 0.0149, but the halfwidth for the same estimate using independent random numbers is larger by a factor of 2.5 at 0.0365. Dierences in Common Independent Output Measures Random Numbers Random Numbers P {T > t0 } 0.0330 +/- 0.0149 0.0330 +/- 0.0365 E[T ] 10.9597 +/- 1.3985 13.1036 +/- 3.8276 P {T1 = T l T < t0 } 0.0370 +/- 0.0124 0.0590 +/- 0.0375 P {T2 = T l T < t0 } 0.0590 +/- 0.0159 0.0640 +/- 0.0381 P {max{T3 , T4 } = T l T < t0 } 0.0960 +/- 0.0197 0.1230 +/- 0.0424 10
Our point estimate of P{T > t0 } is 0.0330 and the sampling error is 0.0149 when using common random numbers. The point estimate indicates that the PM option appears to have higher reliability than the base option, and the fact that 0 is not included in the range 0.0330 +/- 0.0149 = (0.0181, 0.0479) indicates that this dierence is statistically signi"cant at a con"dence level of 95%. Note that our point estimate of P{T > t0 } in Table 3 is identical to the dierence of the point estimates in Table 2. However, the key dierence is that the sampling error of 0.0149 under common random numbers is signi"cantly smaller than the corresponding sampling errors reported in Table 3 for independent random numbers, and signi"cantly smaller than the sampling error for the corresponding absolute performance measures in Table 2. Without hypothesizing a priori whether the base or PM option leads to higher reliability, the question of statistical signi"cance hinges on whether the 95% con"dence interval for P{T > t0 } includes 0. If it does not, the result is statistically signi"cant with the sign of the point estimate of P{T > t0 } determining whether the base or PM option leads to a more reliable system. In this case, a positive dierence indicates the PM option is preferred. Step 6: One-Way Sensitivity Analysis: Sensitivity Plots and Tornado Diagrams From steps 1-4, we have speci"ed a model, restricted attention to key model outputs and inputs, and speci"ed nominal values and ranges for the model inputs. From step 5, we have point estimates, and estimates of sampling error, associated with the models outputs under the models nominal input parameters. Sensitivity plots restrict attention to one or two model outputs at a time and consider a single input parameter. Tornado diagrams restrict attention to one model output at a time and consider multiple inputs. We follow Clemen and Reilly [3] in referring to sensitivity plots and tornado diagrams as one-way sensitivity analysis because we vary one input parameter at a time, holding all other inputs at their nominal values. Figure 2 is a sensitivity plot, showing how system reliability P{T > t0 } changes for the base-option model (without PM) as t0 varies. As t0 grows the system reliability drops. The "gure depicts point estimates along with 95% con"dence intervals on P{T > t0 }. Panel (a) of the "gure shows the results when using common random numbers, and panel (b) shows the same results when using independent random numbers for estimating system reliability. The importance of using common random numbers is evident from the smoothness of the results in panel (a) versus the lack of smoothness in panel (b). Figure 3 is a similar sensitivity plot but for P{T > t0 }, where positive values indicate that the PM option is preferable. The plot indicates a notion of dominance. That is, when other parameters are held at their nominal values, system reliability for the PM option exceeds that of the base option for all t0 values of interest. Figure 4 again shows a sensitivity plot for P{T > t0 }, but now as a function of the reduction in the failure rate k. Here, we see that as k grows the PM option becomes increasingly preferable. We know that at k = 1 the base option should be preferred, but from our simulation results using a sample size of N = 1000, we cannot make this conclusion with statistical signi"cance, as the "gure 11
indicates. (This would change under a larger sample size.) Figures 5 and 6 concern attribution. Here, we suppress the con"dence intervals to avoid clutter, but Tables 2 and 3 provide a sense of the respective 95% con"dence interval halfwidths under common random numbers. These two "gures quantify how the attribution probability of the parallel subsystem 3-4 drops as k grows. Note that sampling error accounts for the dierences in attribution to components 1 and 2 because these components play identical roles and have identical failure rates. 12
Sensitivity Plot: Common Random Numbers 0.90 0.85 Mean 0.80 P[T > t0] 0.75 0.70 0.65 0.60 12 13 14 15 16 17 18 19 20 21 22 23 24 t0 (months) (a) Common Random Numbers Sensitivity Plot: Independent Random Numbers 0.90 0.85 Mean 0.80 P[T > t0] 0.75 0.70 0.65 0.60 12 13 14 15 16 17 18 19 20 21 22 23 24 t0 (months) (b) Independent Random Numbers Figure 2: The y-axis in the "gure is system reliability, P{T > t0 }, and the x-axis is the value of t0 in months for the base-option model. The "gure depicts quantitatively how system reliability drops as the required lifetime of the system grows from its lower bound to its upper bound. The nominal value of t0 is 18 months, and its lower and upper bounds are 12 and 24 months, respectively, as indicated in the "gure. In addition to point estimates of P{T > t0 }, 95% con"dence intervals at each value of t0 are displayed. Panel (a) versus panel (b) of the "gure distinguishes the results when using common random numbers (recommended) versus the more nave approach of using independent random numbers at dierent values of t0 . All estimates are based on a sample size of N = 1000. 13
Sensitivity Plot: Differences - Reliability (to) 0.08 0.07 Mean 0.06 0.05 P[T > t0] 0.04 0.03 0.02 0.01 0.00
-0.01 12 13 14 15 16 17 18 19 20 21 22 23 24 t0 (months)
Figure 3: The "gure is to be read in the same manner as Figure 2 except that P{T > t0 } replaces P{T > t0 } on the y-axis. Sensitivity Plot: Differences - Reliability (k) 0.08 0.07 Mean 0.06 0.05 0.04 P[T > t0] 0.03 0.02 0.01 0.00
-0.01
-0.02 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 k
Figure 4: The "gure is to be read in the same manner as Figure 2 except that P{T > t0 } replaces P{T > t0 } on the y-axis, and k replaces t0 on the x-axis. 14
Sensitivity Plot: PM Option - Attribution (k) 0.60 PM Att 1 0.55 PM Att 2 0.50 P[Ti = T l T < t0] PM Att 3/4 0.45 0.40 0.35 0.30 0.25 0.20 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 k Figure 5: The "gure is to be read in the same manner Figure 2 except that we are plotting the attribution probabilities such as P{T1 = T l T < t0 }. Note that at any vertical line drawn through the three series, the sum of the attribution probabilities is one. Sensitivity Plot: Differences - Attribution (k) 0.20 0.15 0.10 P[Ti = T l T < t0] 0.05 0.00
-0.05
-0.10 Att 1
-0.15 Att 2
-0.20 Att 3/4
-0.25
-0.30 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 k
Figure 6: The "gure is to be read in the same manner Figure 5 except that dierences in attribution probabilities such as P{T1 = T l T < t0 } replace P{T1 = T l T < t0 } on the y-axis. In this case, at any vertical line drawn through the three series the sum of the dierences in the attribution probabilities is zero. 15
A tornado diagram compares the eect of continuously decreasing each input parameter from its nominal value down to its lower bound and increasing the parameter up to its upper bound, and seeing the eect on the models output. Often output measures change monotonically with respect to the inputs. For example, with all other input parameters held constant, increasing the mean time to failure of component 1, 11 , will increase system reliability (P{T > t0 }), increase expected system lifetime (E[T ]), and decrease the probability that a failure is attributable to component 1 (P{T1 = T l T < t0 }). Decreases in 1 1 will have the opposite eect. However, in other cases monotonicity of output is not ensured, and hence we should exercise care that we obtain correct minimum and maximum values of the output as we vary an input over its range. On the x-axis of a tornado diagram we plot the output measure of interest, in this case, system reliability, P{T > t0 }, expected time to system failure, E[T ], or their dierences, P{T > t0 } and E[T ]. The y-axis stacks bars with the range of these outputs for each input parameter of interest. The output under the nominal values of the inputs is highlighted. The horizontal bars for the input parameters are ordered by sensitivity, with the longest bar, i.e., most sensitive input parameters on top. Note the importance of having selected commensurate ranges for the input variables in step 4, as these now aect which parameters are seen as most important. Panels (a) and (b) of Figure 7 display tornado diagrams for the base-option model for P{T > t0 } and E[T ]. Panels (a) and (b) of Figure 8 are analogous but for the PM option. Panels (a) and (b) of Figure 9 display P{T > t0 } and E[T ], where positive values indicate the PM option has higher reliability and longer expected system lifetime. Panel (a) of Figure 7 indicates that system reliability under the base option is most sensitive to the value of t0 , followed by the mean failure times of components 3 and 4 and then the mean times for components 1 and 2. System reliability under the base option is not aected by input parameters k and t. For E[T ], panel (b) of Figure 7 indicates the most sensitive parameters are the mean failure times of components 3 and 4 followed by the same times for components 1 and
- 2. Output E[T ] is not aected by input parameters k, t, or t0 . The results for the PM option in Figure 8 are similar except that we see the importance of parameter k. Interestingly, the results are relatively insensitive to the duration of the PM interval, t, although this changes slightly in Figure 9 when examining the results for P{T > t0 }. For this reason, we do not show sensitivity plots with respect to t here, although we revisit sensitivity to t in step 8 (in spider plots) and step 9 (in a two-way sensitivity analysis with k).
Sampling errors are more easily displayed on sensitivity plots than on tornado diagrams. Still, in Figures 7-9 we display horizontal error bars at the two extremes. Sensitivity plots also have the advantage that multiple outputs can be plotted simultaneously. That said, tornado diagrams can display many input parameters simultaneously and are widely used to assess to which input parameters an output is most sensitive. 16
(a) P{T > t0 } (b) E[T ] Figure 7: The "gure depicts tornado diagrams for P{T > t0 } and E[T ] for the base-option model. Point estimates along with 95% con"dence intervals are displayed. The "gure indicates that system reliability is most sensitive to the value of t0 , followed by the mean failure times of 1 1 3 and 4 and then by the mean times 1 and 2 . Expected system lifetime is most sensitive to 3 and 1 1 1 1 4 and then 1 1 1 and 2 . Note that the color shading indicates whether a high or low value of the input parameter corresponds to the change in the output. Thus, higher values of t0 lead to lower values of system reliability, but higher component mean lifetimes lead to higher system reliability. Again, estimates and error bars are calculated based on a sample size of N = 1000. 17
Tornado Diagram: PM Option - P[T > t0] 0.66 0.68 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84 0.86 0.88 t0 k 1/3 1/2 1/4 1/1 t Low High (a) P{T > t0 } Tornado Diagram: PM Option - E[T] 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80 k 1/3 1/4 1/2 1/1 t t0 Low High (b) E[T ] Figure 8: The "gure depicts tornado diagrams for P{T > t0 } and E[T ] for the PM option, and is to be read in the same manner as Figure 7. 18
Tornado Diagram: Difference - P[T > t0]
-0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 k
t0 1/4 1/3 t 1/1 1/2 Low High (a) P{T > t0 } Tornado Diagram: Difference - E[T]
-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 k
1/3 1/2 1/4 1/1 t t0 Low High (b) E[T ] Figure 9: The "gure depicts tornado diagrams for P{T > t0 } and E[T ] where positive values favor the PM option. 19
Step 7: Uncertainty Quanti"cation Plots An important part of uncertainty quanti"cation (UQ), and a part that distinguishes it from routine sensitivity analysis, concerns propagating a probability distribution placed on one or more input parameters through the nonlinear function represented by a simulation model and characterizing the resulting probability distribution on an output measure. We emphasize that the probability distribution we speak of here is a probability model that we place on input parameters and not the Monte Carlo sampling-based error we reference above. (That said, as elsewhere, we also capture sampling-based error here, too.) We call a graphical plot of the resulting probability distribution a UQ plot, regardless of whether it is expressed as a cumulative distribution function (cdf) or a probability density function (pdf). This idea is closely related to the sensitivity plots we form in step 6, except that we now embed information associated with the probability distribution placed on the input parameters. We begin by focusing on the case when a univariate distribution is placed on a single input parameter, and we then turn to UQ plots when multivariate distributions are placed on input parameters. When constructing a sensitivity plot, the y-axis is the output parameter, and the x-axis is the input parameter. For a sensitivity plot we typically form a uniform grid over the range of the input parameter values, e.g., over the input ranges that Table 1 speci"es. A cdf-based UQ plot is a plot of the cumulative distribution function of the output measure. We also form pdf-based UQ plots. In both cases we form estimates of these function based on sampling, where the sampling is done in a manner we make precise below. For a cdf UQ plot, the x-axis contains levels of the output measure, the y-axis contains probabilities, and the probability distribution on the input parameter is implicitly encoded in the result. We again use our example to make this idea concrete. Suppose that the improvement factor, k, has a continuous uniform random variable on the interval (1, 5) speci"ed in Table 1. Figure 10 contains UQ plots of E[T l k] and P{T > t0 l k} for the PM-option model. The two panels of the "gure contain estimates of the cdf-based UQ plots for both these two outputs. As k grows the probability that a system failure is due to a failure of the parallel subsystem of components 3-4 drops. As a result, we see both cdfs grow quickly towards one for large values of E[T l k] and P{T > t0 l k} because there is a large probability mass for k associated with little improvement in these values. Figure 11 is similar, except that we now show cdfs for E[T l k] and P{T > t0 l k} rather than E[T l k] and P{T > t0 l k}. Finally, Figure 12 shows the pdfs for E[T l k] and E[T l k]. We do not show analogous pdfs for P{T > t0 l k} and P{T > t0 l k} because the estimates have excessive sampling error. Developing good pdf estimates for P{T > t0 l k} and P{T > t0 l k} would require a larger sample size. 20
UQ Plot: PM Option - CDF E[Tlk] 1.00 0.90 Mean 0.80 UL 0.70 CDF of E[Tlk] LL 0.60 0.50 0.40 0.30 0.20 0.10 0.00 40 45 50 55 60 65 70 75 80 E[Tlk] (a) E[T l k] cdf UQ Plot: PM Option - P[T > t0lk] 1.00 0.90 Mean 0.80 UL CDF of P[T > t0 lk] 0.70 LL 0.60 0.50 0.40 0.30 0.20 0.10 0.00 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83 0.84 P[T > t0lk] (b) P{T > t0 l k} cdf Figure 10: The "gure depicts UQ plots which consist of estimates of the cdf of the corresponding output measures when the improvement factor, k, is a uniform random variable on the interval (1, 5). Point estimates as well as a 95% con"dence envelope are plotted. As k grows the probability that a system failure is due to the parallel subsystem of components 3-4 shrinks. As a result, we see both cdfs grow quickly towards one for large values of E[T l k] and P{T > t0 l k}. 21
UQ Plot: Difference - CDF E[Tlk] 1.00 0.90 Mean 0.80 UL 0.70 CDF of E[Tlk] LL 0.60 0.50 0.40 0.30 0.20 0.10 0.00
-5 0 5 10 15 20 25 30 35 E[Tlk]
(a) E[T l k] cdf UQ Plot: Difference - P[T > t0lk] 1.00 0.90 Mean 0.80 CDF of P[T > t0lk] UL 0.70 LL 0.60 0.50 0.40 0.30 0.20 0.10 0.00
-0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 P[T > t0lk]
(b) P{T > t0 l k} cdf Figure 11: The "gure is to be read as Figure 10 except that we now show cdfs for E[T l k] and P{T > t0 l k} rather than E[T l k] and P{T > t0 l k}. 22
UQ Plot: PM Option - PDF E[Tlk] 0.08 Mean 0.07 UL LL 0.06 PDF of E[Tlk] 0.05 0.04 0.03 0.02 0.01 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 E[Tlk] (a) E[T l k] pdf UQ Plot: Difference - PDF E[Tlk] 0.08 Mean 0.07 UL LL 0.06 PDF of E[Tlk] 0.05 0.04 0.03 0.02 0.01
-5 0 5 10 15 20 25 30 E[Tlk]
(b) E[T l k] pdf Figure 12: The "gure depicts UQ plots which consist of estimates of the pdf of E[T l k] and E[T l k] when the improvement factor, k, is a uniform random variable on the interval (1, 5). 23
In our example, when forming UQ plots of E[T l k] and P{T > t0 l k} we regard the other six input parameters, 1 , . . . , 4 , t, and t0 , as deterministic parameters, and the failure times of the four components, T1 , . . . , T4 as random variables. Our sampling consists of drawing N = 1000 independent and identically distributed (i.i.d.) observations of the four-tuple (T1 , . . . , T4 ). In this one-dimensional setting we form the 1%, 2%, 3%, . . . , 99% percentiles of the distribution of k, using its distribution, and we then use our N = 1000 i.i.d. observations of (T1 , . . . , T4 ) to estimate E[T l k = k ], for = 0.01, 0.02, . . . , 0.99, where k denotes these percentiles. Although we describe conditioning on evenly-spaced quantilesevenly spaced in terms of probabilityit may be desirable to have a "ner grid in regions where the function changes most rapidly. Again, we emphasize the importance of using common random numbers in forming UQ plots such as those in Figures 10-12. Obvious alternatives to what we have just sketched are also possible, appropriate, and even necessary. (It also important to recognize what is inappropriate and we point to that, too, below.) For example, we could regard the other six parameters as random variables instead of "xing them at their nominal values and sample them in the same way we sample (T1 , . . . , T4 ), while still conditioning on k = k to form estimates of E[T l k = k ]. Clarity in exposition should indicate what precisely the expected-value operator is averaging over. In another alternative, we could also sample from ks distribution instead of conditioning on its quantiles in order to form a UQ plot. When sampling k, it is important to distinguish this sampling from that for (T1 , . . . , T4 ). Speci"cally, we could use one sample size Nuq for k and form k i , i = 1, . . . , Nuq from ks distribution. For each of these samples we then compute, or rather estimate E[T l k = k i ], i = 1, . . . , Nuq , where each estimate averages over the N sampled realizations of (T1 , . . . , T4 ). We then use the estimates of E[T l k = k i ], i = 1, . . . , Nuq , to form the types of plots in Figure 10. We did not use this sampling-based method in forming the UQ plots of Figure 10 because it is more ecient in the one-dimensional setting to condition on the quantiles as we describe above. However, this sampling-based approach is necessary to form a UQ plot when bivariate, or higher dimensional multivariate, distributions are placed on input parameters. For example, if we place a bivariate distribution on the PM time-reduction factor pair, (t, k), then such a bivariate distribution has no notion of quantiles, and so the one-dimensional procedure does not have a bivariate analog. We must sample. Importantly, we do not have to sample (T1 , . . . , T4 ) from its underlying distribution. If we have a variance reduction scheme that forms, e.g., unbiased estimates of E[T l k = k i ] then that sampling scheme can be used. However, it is important to recognize that we must sample from the true underlying distribution of (t, k). If we use a distribution-altering variance reduction scheme to sample from (t, k) when averaging out those parameters, that scheme cannot be used when forming a UQ plot. Such altered schemes are designed to reduce variance and hence would yield misleading plots, indicating, e.g., a pdf that is too narrow about the mean value of the output measure. 24
Step 8: One-Way Sensitivity Analysis: Spider Plots A spider plot is an x-y graph, in which the x-axis depicts changes in the input parameters and the y-axis captures corresponding changes in the models output measure. Like a tornado diagram, a spider plot involves multiple input parameters and a single output variable. The output variable is typically expressed in its natural units. For example, we express changes in E[T ] in months and we express changes in P{T > t0 } as a unitless value between 0 and 1, where the changes are relative to estimates under the nominal value of the parameters. In order to allow the x-axis to simultaneously represent multiple input parameters, which are on dierent scales with dierent units, there are two possibilities. One possibility is to express percentage changes in the input parameters from their nominal values. The second possibility is to express changes as multiples of the standard deviation of the input parameters, when those input parameters are governed by probability distributions. In either case, the magnitude we vary the parameters is determined by the reasonable and commensurate ranges we have speci"ed in step 4 of the analysis, e.g., in Table 1. In the former case, if the nominal value of the input parameter is zero, then a second x-axis must be added. A tornado diagram can include a larger number of input variables than a spider plot. A spider plot allows displaying about seven input parameters before it becomes cluttered. If the output variable is monotonic (increasing or decreasing) in an input parameter then we only need to estimate the models output at the lower bound, nominal value, and upper bound of the input parameter. A spider plot requires estimating the models output at enough values of each input parameter that a seemingly continuous plot of (x, y) pairs can be formed. A spider plot contains more information than a tornado diagram. The tornado diagrams endpoints denote the endpoints of the spider plot, but the spider plot also speci"es changes in the output at intermediate values, as does a sensitivity plot. Again like a sensitivity plot, we can assess whether changes in the output are linear or nonlinear with respect to changes in the input. Spider plots can contain point estimates or 95%, say, con"dence intervals on those changes. (In the latter case, we may need to reduce the number of input parameters simultaneously displayed.) Figure 13 displays two spider plots for our example for E[T ] and P{T > t0 } (y-axis) for the base-option model as a function of percentage changes in the input parameters (x-axis). Panel (a) shows a spider plot for E[T ] while panel (b) shows the spider plot for P{T > t0 }. Figure 14 displays spider plots for the PM option for our example, and Figure 15 displays spider plots for E[T ] and P{T > t0 }. Note that we choose not to include con"dence limits for each parameter displayed in the graph because of the clutter they induce. Also note that the range of the x-axis is determined by the nominal values and the lower and upper bounds speci"ed in Table 1. Importantly, the input parameters are not each varied the same percentage. Rather, the limits are those speci"ed in Table 1. Qualitatively, the "gures are similar to the tornado diagrams as to the most sensitive input parameters. However, Figures 13-15 are more insightful as to the rate of change and any associated nonlinearities in the change. 25
Spider Plot: Base Option - E[T] 6 4 2 1/1 Change in E[T] 1/2 0 1/3
-2 1/4
-4
-6
-8
-50 -40 -30 -20 -10 0 10 20 30 40 50
% Change in Input Parameters (a) E[T ]
Spider Plot: Base Option - P[T > t0] 0.10 0.05 Change in P[T > t0] 1/1 0.00 1/2 1/3
-0.05 1/4 t0
-0.10
-0.15
-50 -40 -30 -20 -10 0 10 20 30 40 50
% Change in Input Parameters (b) P{T > t0 }
Figure 13: The "gure depicts two spider plots for the base option in our example. The plot in panel (a) shows E[T ] (y-axis) as a function of percentage changes in the input parameters (x-axis). The plot in panel (b) is identical to the one in panel (a) except that the y-axis is P{T > t0 }. 26
Spider Plot: PM Option - E[T] 20 15 1/1 10 Change in E[T] 1/2 5 1/3 0 1/4
-5 k
-10 t
-15
-50 -25 0 25 50 75 100 125 150
% Change in Input Parameters (a) E[T ]
Spider Plot: PM Option - P[T > t0] 0.08 0.06 1/1 0.04 Change in P[T > t0] 1/2 0.02 1/3 0.00 1/4
-0.02 k
-0.04
-0.06 t
-0.08 t0
-0.10
-50 -25 0 25 50 75 100 125 150
% Change in Input Parameters (b) P{T > t0 }
Figure 14: The "gure reads as Figure 13 except that it is for the PM option in our example. 27
Spider Plot: Difference - E[T] 20 15 1/1 10 Change in E[T] 1/2 5 1/3 0 1/4
-5 k t
-10
-15
-50 -25 0 25 50 75 100 125 150
% Change in Input Parameters (a) E[T ]
Spider Plot: Difference - P[T > t0] 0.03 0.02 1/1 Change in P[T > t0] 0.01 1/2 1/3 0.00 1/4
-0.01 k
-0.02 t
-0.03 t0
-0.04
-50 -25 0 25 50 75 100 125 150
% Change in Input Parameters (b) P{T > t0 }
Figure 15: The "gure reads as Figure 13 except that it is for the dierences in performance for our example. 28
Step 9: Two-way Sensitivity Analysis Two-way sensitivity graphs allow for visualizing the interaction of two or more input variables. While such analysis can be more dicult to perform, it can provide valuable insight. For our example, Figure 16 depicts the eect of simultaneous changes in the duration of PM (t) and the factor by which the PM reduces the failure rate of component 3 (k) on P{T > t0 }. Panel (a) contains only the point estimate, and panel (b) contains both the point estimate and the con"dence limits, allowing us to see the indierence zone where neither option is statistically better than the other. Figure 17 is similar to Figure 16, except the parameters of interest are t0 and k. Two-way sensitivity analyses can be extended to include more than two decision alternatives, and in such cases the plots typically partition the space into three or more regions in which each alternative is preferred. It is also possible to form a three-dimensional plot of an output variable (e.g., the dierence in system reliability P{T > t0 }) as a function of two input variables (e.g., t and k), although we do not pursue that here. 29
Two-Way Sensitivity Plot: Difference - P[T > t0] 5 Base Option 4 Preferred 3 t 2 PM Option Mean Preferred 1 0 1 2 3 4 5 k (a) Point Estimate Two-Way Sensitivity Plot: Difference - P[T > t0] 5 Base Option 4 Preferred 3 t Mean 2 PM Option LL Preferred UL 1 0 1 2 3 4 5 k (b) Point Estimate and Con"dence Limits Figure 16: The "gure depicts a two-way sensitivity plot for the input parameters governing the duration of PM, t, and the factor by which the failure rate of component 3 is reduced, k. When t is small and k is large, the PM option is preferred, and when t is large and k is small the base option is preferred. The "gure quanti"es this notion with the two-dimensional analog of a threshold analysis. 30
Two-Way Sensitivity Plot: Difference - P[T > t0] 24 PM Option 22 Preferred 20 t0 18 16 Base Option Mean Preferred 14 12 1 2 3 4 5 k (a) Point Estimate Two-Way Sensitivity Plot: Difference - P[T > t0] 24 PM Option 22 Preferred 20 t0 18 Mean UL 16 Base Option LL 14 Preferred 12 1 2 3 4 5 k (b) Point Estimate and Con"dence Limits Figure 17: The "gure depicts a two-way sensitivity plot for the input parameters governing the minimum time the system must be operational, t0 , and the factor by which the failure rate of component 3 is reduced, k. When t0 is small and k is small, the base option is preferred, and when t0 is large and k is large the PM option is preferred. 31
Step 10: Metamodels & Design of Experiments The possible pitfalls of changing only one input parameter at a time are well documented in the literature (see, e.g., [8, 10]) and include the fact that interaction eects between input parameters are lost. Graphical sensitivity analyses become more dicult when moving past a one- or two-dimensional analysis, but we can form a metamodel (which is also called a response surface, an emulator, or a surrogate model) and carry out an experimental design to "t that metamodel. Recall from step 1 that we use y to denote the simulation models output and we use x = (x1 , . . . , xn ) to denote the simulation models input. For simplicity we focus on a single output measure. Among the simplest metamodels typically postulated are polynomial regression models of low degree; e.g.,
n n n y = 0 + k xk + k,k xk xk + . (3) k=1 k=1 k =k+1 To "t the parameters 0 , 1 , . . . , n , and 1,2 , . . . , 1,n , . . . , n1,n , we use an experimental design to specify M , say, input parameter vectors (xi1 , . . . , xin ), i = 1, . . . , M , coupled with the corresponding estimated simulation output y i . The corresponding error terms, i , are assumed to be independent and normally distributed with mean zero, at least in the simplest approach. So far, we have em-phasized the importance of using common random numbers in carrying out our analyses. However, the assumption that the error terms i are independent requires that we draw independent Monte Carlo samples at each design point. When estimating a dierence, we still use common random numbers for the base and PM options at that design point.
We could add higher-order cross terms to equation (3). This can improve the quality of the "t, but we need to take care that we have an adequate number of observations relative to the larger number of model parameters so the model is not over-"t. Perhaps more importantly, a thorough understanding of a higher-degree polynomial regression model can be challenging. Table 1 speci"es seven input parameters for our example, but only "ve of those input parameters matter for the base option. (The PM parameters of t and k are not relevant for the base option.) In this case, we can seek to explain system reliability y = P{T > t0 } as a relatively simple function of 1 1 1 1 1 , 2 , 3 , 4 , and t0 . A full factorial design with two levels in this case would involve 25 = 32 combinations of these values placed at their lower and upper bounds from Table 1. When we also include the nominal case for the values of the parameters we have 35 = 243 combinations. If we instead have seven input parameters for the PM option, these values become 27 = 128 and 37 = 2187. Note that the preceding values represent the number of design points in a speci"c experimental design. In order to estimate higher-order eects and interaction terms, multiple replications are typically run at each design point. So for our example, at each design point, we could estimate y = P{T > t0 } using a sample size of N = 1000, and we could replicate that point estimate Nr = 3 times. For our example, we present the results of a regression metamodel, where the performance measure of interest is the dierence in reliability, P{T > t0 }, where again, positive values favor 32
the PM option. We use a 37 full factorial design, where the three levels are the minimum, nominal, and maximum values in Table 1. At each design point, we perform Nr = 3 independent replications of the simulation using a sample size of N = 1000 for each replicate. Table 4: Experimental design parameters for the replicated 37 factorial design and overall model "t statistics. The F -Statistic measures whether any of the input variables in any combination with one another have a statistically signi"cant eect on the response (P{T > t0 }). The adjusted R2 value assesses the overall goodness-of-"t of the model to the simulation output.
# of Design Points 2,187
# of Replications Per Design Point: 3 Total Sample Size: 6,561 Residual Standard Error: 0.01295 Residual Degrees of Freedom: 6,532 Adjusted R2 Value: 0.8895 F -Statistic: 1,886 Degrees of Freedom 1: 28 Degrees of Freedom 2: 6,532 p-value: <2.20E-16 From the information in Table 4, we can see that based on the adjusted R2 value of 0.8895, we have a good (not excellent) "t to the simulation output. An excellent "t is categorized by an adjusted R2 value of at least 0.90. As we indicate above, in a 37 full factorial experimental design, there are 2,187 design points; i.e., 2,187 unique combinations of input parameters. We replicate these design points three times, for a total size of 6,561. Each of these 6,561 estimates of P{T > t0 } is based on a sample size of N = 1000. The F -Statistic and p-value are used to test the hypothesis that all of the coecients k and k,k , with the exception of the intercept 0 , are zero.
If this is the case, none of the input parameters (i.e., factors in the terminology of regression and experimental design) or interactions among these parameters are signi"cant in estimating the response variable, P{T > t0 }. The p-value measures the statistical signi"cance of the result. For the F -Test, we see that the p-value is less than 2.2 x 1016 , which is much less than the standard signi"cance level of 0.05. So for this F -Test we can conclude that at least one coecient in our regression model is signi"cantly dierent from zero, and hence is a signi"cant factor in estimating P{T > t0 }. This, coupled with a reasonably good adjusted R2 value of 0.8895 suggests that we have a relatively good "t to the simulation output and the "t parameters, other than the intercept, play a signi"cant role in estimating P{T > t0 }. The next step in the analysis is to examine which regression coecients, k and k,k , are statis-tically signi"cant in our model. This is a notion of attribution in that we are attempting to attribute a change in the response variable to changes in the input variable and characterize the strength of the relationship. Table 5 presents the estimates for the coecients in the regression model (3) as well as their standard errors, and the associated p-values for determining the signi"cance levels of 33
individual parameters. Table 5: Results of "tting linear regression model (3) for P{T > t0 }. We exclude all terms higher than second order terms as equation (3) indicates. Using the p-value we can determine which eects and two-way interactions have a signi"cant eect in estimating P{T > t0 }. Coecient Standard Signi"cant Coecients Parameters t-value p-value Estimates Error at = 0.05 0 Intercept -4.374E-02 8.909E-03 -4.910 9.34E-07 Yes 1 1 1 -4.420E-06 3.132E-05 -0.141 8.88E-01 No 2 1 2 -4.166E-05 3.132E-05 -1.330 1.84E-01 No 3 1 3 2.448E-04 8.680E-05 2.820 4.82E-03 Yes 4 1 4 4.089E-04 8.680E-05 4.711 2.52E-06 Yes 5 k 1.480E-02 9.005E-04 16.434 2.00E-16 Yes 6 t -2.160E-02 1.801E-03 -11.993 2.00E-16 Yes 7 t0 2.959E-03 2.816E-04 10.507 2.00E-16 Yes 1,2 1 1
- 2 1
7.682E-09 9.591E-08 0.080 9.36E-01 No 1,3 1
- 1 1
3 -3.916E-07 2.398E-07 -1.633 1.02E-01 No 1,4 1 1
- 1 4 -3.837E-07 2.398E-07 -1.600 1.10E-01 No 1
1,5 1
- k 1.009E-05 2.398E-06 4.206 2.63E-05 Yes 1,6 1 1
- t -6.564E-06 4.796E-06 -1.369 1.71E-01 No 1,7 11
- t0 1.901E-06 7.993E-07 2.378 1.74E-02 Yes 2,3 2
- 1 1
3 -4.407E-07 2.398E-07 -1.838 6.61E-02 No 2,4 1 2
- 1 4 -1.680E-08 2.398E-07 -0.070 9.44E-01 No 1
2,5 2
- k 1.061E-05 2.398E-06 4.427 9.73E-06 Yes 2,6 1 2
- t -5.576E-06 4.796E-06 -1.163 2.45E-01 No 2,7 12
- t0 3.070E-06 7.993E-07 3.842 1.23E-04 Yes 3,4 3
- 1 1
4 1.319E-05 5.994E-07 21.998 2.00E-16 Yes 3,5 1 3 *k -2.303E-04 5.994E-06 -38.416 2.00E-16 Yes 3,6 1 3
- t -8.428E-05 1.199E-05 -7.030 2.28E-12 Yes 3,7 13
- t0 -3.768E-05 1.998E-06 -18.858 2.00E-16 Yes 4,5 1 4 *k -2.378E-04 5.994E-06 -39.675 2.00E-16 Yes 4,6 1 4
- t 2.797E-04 1.199E-05 23.329 2.00E-16 Yes 4,7 14
- t0 -5.333E-05 1.998E-06 -26.691 2.00E-16 Yes 5,6 k
- t -1.725E-03 1.199E-04 -14.390 2.00E-16 Yes 5,7 k
- t0 1.183E-03 1.998E-05 59.200 2.00E-16 Yes 6,7 t
- t0 3.536E-04 3.996E-05 8.847 2.00E-16 Yes Table 5 provides several insights, and can be viewed as typical regression output. In the "rst two columns of the table, we list all of the coecients in the model, their corresponding input parameters, and interactions among the input parameters. For example, 0 is the coecient representing the intercept, 1 represents the coecient of 1 1 , and 5,7 represents the interaction between k and t0 . The third column presents point estimates of each of the regression coecients and the fourth column presents the standard error associated with these coecients. The ratio of 34
the estimate to the standard error is the t-value for the t-test, which tests the hypothesis that each individual parameter takes value zero so that rejecting the null hypothesis indicates the parameter is signi"cant. The p-values for the t-tests for each of these terms (column 6) can be interpreted in same manner as the p-value for the F -test we describe above, and can be used to determine the statistical signi"cance of each term in the regression model individually, whereas the F -test determines the signi"cance of the parameters as a group. The positive signs on the coecients for k and t0 are consistent with the sensitivity plots in Figures 3 and 4 and the spider plot in Figure 15 and, of course, with intuition. Similarly, the negative coecient for t is consistent with intuition and Figure 15: As the time required to carryout PM grows, the bene"t of the PM option drops. The signs and magnitude of other coecients are more subtle. The positive coecient for 1 3 is counterintuitive: As the reliability of component 3 grows, the bene"t from PM should shrink, not grow. However, examining the sign and relative magnitude of the coecients for 1 1 3 and 3
- kand knowing the nominal value of k is 2we see that as 1
3 grows, the regression estimate of P{T > t0 } indeed shrinks. This holds except for values of k 1, on the boundary of the experimental design, where the quality of the regression "t is likely poorer. We can see from the information in Table 5, that several of the two-factor interactions terms are signi"cant, including all interaction terms involving k and t0 , indicating that these parameters have a signi"cant eect on the response variable P{T > t0 }. This is again consistent with what we learned from the tornado diagram in Figure 9 and the spider plot of Figure 15. In addition, we see that the interaction term k
- t0 has the greatest t-statistic among all terms, and we can conclude that it is one of the most signi"cant contributors to the estimate of P{T > t0 }. The fact that this nonlinearity is important, and that the sign of k
- t0 s coecient is positive, is not surprising given the two-way sensitivity plot in Figure 17 for k and t0 .
To illustrate the value of this regression output, we demonstrate how to use this regression equation to predict the value of P{T > t0 }, without the need to rerun the simulation model. This is especially important for large-scale stochastic simulation models for which it may take several days or even weeks to run a designed experiment and collect the type of output we collected for this illustrative example. If the statistical experiment is planned well, the regression metamodel can be used to predict performance measures without the need for rerunning the simulation model. However, we note that using this approach, we should not use values of the input parameters outside the bounds of those used in the experimental design, in our case the minimum parameter values from Table 1. And, as we have already seen, near the boundary of the design, the regression "t may degrade. In addition, we should examine the goodness-of-"t statistics we describe above to ensure the statistical model is adequate before relying on its predictive power. Table 6 applies our regression metamodel to the nominal values for all parameters, and forms an estimate of P{T > t0 } without the need to run the simulation model. Table 7 presents the results of simply using the simulation model to estimate P{T > t0 } and its associated error, and we can compare these estimates to the value obtained by using the regression metamodel to determine the validity 35
of metamodel. Table 6: An example of applying the results of the regression model to the nominal values of the input parameters in order to predict P{T > t0 }. Note that using this method, we no longer need to run the simulation model to predict P{T > t0 } as a function of any combination of input parameters . The adjusted R2 value from Table 4 is 0.8895, and in practice this value is acceptable, although 0.90 is typically considered the threshold for a excellent "t. Coecient Parameter Coecients Parameters Contribution Estimates Values 0 -4.374E-02 Intercept None -0.0437 1 -4.420E-06 1 1 200 None 2 -4.166E-05 1 2 200 None 3 2.448E-04 1 3 50 0.0122 4 4.089E-04 1 4 50 0.0204 5 1.480E-02 k 2 0.0296 6 -2.160E-02 t 1 -0.0216 7 2.959E-03 t0 18 0.0533 1,2 7.682E-09 1 1
- 2 1
40000 None 1,3 -3.916E-07 1 1
- 3 1
10000 None 1,4 -3.837E-07 1
- 1 1
4 10000 None 1,5 1.009E-05 1 1
- k 400 0.0040 1,6 -6.564E-06 1 1
- t 200 None 1,7 1.901E-06 11
- t0 3600 0.0068 2,3 -4.407E-07 1 2
- 3 1
10000 None 2,4 -1.680E-08 2
- 1 1
4 10000 None 2,5 1.061E-05 1 2
- k 400 0.0042 2,6 -5.576E-06 1 2
- t 200 None 2,7 3.070E-06 12
- t0 3600 0.0111 3,4 1.319E-05 1 3
- 4 1
2500 0.0330 3,5 -2.303E-04 1 3 *k 100 -0.0230 3,6 -8.428E-05 1 3
- t 50 -0.0042 3,7 -3.768E-05 13
- t0 900 -0.0339 4,5 -2.378E-04 1 4 *k 100 -0.0238 1
4,6 2.797E-04 4
- t 50 0.0140 4,7 -5.333E-05 14
- t0 900 -0.0480 5,6 -1.725E-03 k
- t 2 -0.0035 5,7 1.183E-03 k
- t0 36 0.0426 6,7 3.536E-04 t
- t0 18 0.0064 Estimate of P{T > t0 } 0.03591 36
Table 7: Results for P{T > t0 } from running the simulation model with all input parameters at nominal levels with sample size N = 1000. P{T > t0 } 0.0330 95% CI Halfwidth 0.0149 95% CI Lower Limit 0.0181 95% CI Upper Limit 0.0479 In Table 6 we again provide the estimates of the coecients and the parameters and interactions they represent. The fourth column of this table contains the parameter values for each of the associated coecients for the nominal case, with the exception of the intercept which is not directly linked to any parameter or interaction. The "rst seven values in this column match the nominal parameter values given in Table 1. The interaction terms are simply the products of these input parameter values. For example, the nominal value of t = 1, and the nominal value of k = 2, and thus the parameter value for the interaction t*k = 2, as shown in the table. As shown in regression model (3), an estimate of P{T > t0 } can be formed by computing the product of the parameter (or interaction term) values and the coecient, and summing those values along with the value of the intercept. The "nal column in Table 6 is the product of the coecient and parameter value. Summing all the values in this column provides us with an estimate of P{T > t0 } = 0.03591 when all parameters are at their nominal levels. It is useful to compare this estimate with a point estimate and associated 95% con"dence limits of P{T > t0 } using the simulation model. Table 7 presents results from the simulation model, and we see that the estimate of P{T > t0 } is 0.0330 +/- 0.0149 = (0.0181, 0.0479). We can see that our regression metamodel estimate of 0.0359 is within the 95% con"dence limits, and is dierent from the simulation point estimate of 0.0330 by less than 10%. This suggests that our regression metamodel can be used in lieu of running the simulation under appropriate circumstances. In addition to predicting the values of performance measures, we can also use the results of the designed experiment to construct two-way interaction plots that describe how a given response changes as a function of two input parameters. Figures 18-21 show two-way interaction plots for pairs of input parameters, where the response variable is again P{T > t0 }. For example, in Figure 18, we see how P{T > t0 } changes when all input parameters are held constant except for k and 1 3 . First, we see that when k = 1, using the point estimates, we prefer the base option regardless of the value of 1 3 . However, for values of k greater than one (namely 3 and 5), we prefer the PM option regardless of the value of 1 1 3 . When 3 is at its minimum value of 30, changes in the repair factor k have a more signi"cant eect on P{T > t0 } than when 1 3 is at the nominal or maximum levels. Figure 19 shows similar results for k and t0 . The fact that P{T > t0 } grows with k and t0 is as expected, as is the ampli"cation of the eect of growing k for larger values of t0 . Figures 20 and 21 depict analogous results for the respective pairs (k, t) and (t0 , t). 37
Figure 18: The "gure is a two-way interaction plot for k and 1 3 , where the response variable is P{T > t0 }. Figure 19: The "gure is a two-way interaction plot for k and t0 , where the response variable is P{T > t0 }. 38
Figure 20: The "gure is a two-way interaction plot for k and t, where the response variable is P{T > t0 }. Figure 21: The "gure is a two-way interaction plot for t0 and t, where the response variable is P{T > t0 }. 39
Running full factorial designs, as the number of input parameters grows large, and the underly-ing simulation model is computationally expensive, quickly becomes intractable. The way forward is to employ fractional experimental designs and/or to attempt to reduce the dimension of the input vector. Fractional factorial designs use fewer design points, and like the analysis we give above, disregard higher order interaction terms in favor of estimating main eects and two-way interactions. There is a large literature on this topic, and it is not our goal here to review this in detail. See, for example, the survey in [9] and references cited therein. 4 Further Discussion Implicit in much of our discussion in steps 5-9 is the notion of threshold analysis. For our example, the PM option is preferred under the nominal values of the input parameters. Understanding how much an input parameter would need to change in order to change that assessment has been important in our discussions of one- and two-parameter sensitivity plots, tornado diagrams, UQ plots, and spider plots. This notion is even embedded in our de"nition of the output performance measure P{T > t0 } P{TP M > t0 } P{Tbase > t0 }. Strictly speaking, the notion of a threshold value applies when estimates of the output variables are precise. When they contain sampling error, or they are uncertain because of uncertainties in the input parameters, a more nuanced analysis is needed. For example, in a sensitivity plot for P{T > t0 } as a function of k, we examine whether 95% con"dence intervals include zero to understand whether the PM option is preferred, the base option is preferred, or whether, based on sampling error, we cannot assess which is preferable. This third characterization is made if the con"dence interval contains zero. In that case, we are in an indierence zone in which we cannot assess whether the PM or base option is preferred. When we place a prior distribution on k, we can obtain a probability distribution governing the random variables P{T > t0 l k}, and E[T l k], and we can assess the probability that one of these outputs is positive (favoring the PM option) or negative (favoring the base option). In some cases, as we range an input parameter, the preference for one option over an alternative does not change. In this case we have a dominance relationship. As Figure 3 illustrates, we prefer the PM option over the base option for all values of the input parameter t0 of interest. It is important to note that we may have considerable variability in an output measure, but if we have a dominance relationship then this variability is of secondary interest. Of primary interest is whether the decision we would make changes. Perhaps the foremost caveat when performing a one-way sensitivity analysis of a trusted model is to assess whether it makes sense to vary the input parameters one at a time. If two or more input parameters depend on an unstated auxiliary factor, this may not be valid. For example, it may be that components 1 and 2 in our example are identical but we are unsure of the associated failure rate. In this case, we should have 1 = 2 , i.e., we should replace the two input parameters 1 and 2 with a single parameter. In other cases, the dependency between two parameters is not so simple. For example, if 1 and 2 are random variables they may have a dependent joint 40
probability distribution in which the correlation is positive (but not perfect). In this case, we could view that correlation coecient between 1 and 2 as an input parameter to be varied or we could form a UQ plot and to characterize P{T > t0 l (1 , 2 )}. Another caveat in our one-at-a-time sensitivity analysis concerns the notion that all input parameters take their nominal values except for one. This can mask the eect of cross terms. The level of one input variable departing from its nominal value may amplify the eect of changes in another input variable. The purpose of the metamodel analysis in step 10 is to unmask precisely such interactions. Whether done by a formal expert elicitation or an informal scheme, the nominal values of the parameters and their ranges are typically based on expert opinion and hence subject to well-known biases related to anchoring, over con"dence, etc. Particularly relevant for GSI-191 analysis are diculties in assessing rare-event probabilities. For example, see the discussions in Tversky and Kahneman [21] and OHagan et al. [15], and also see Kynn [11]. Model uncertainty is an often neglected part of sensitivity analysis. A simple form of model uncertainty for our example concerns the distributional assumption on the four-tuple of failure times (T1 , . . . , T4 ). We have assumed these four failure times to be independent and to have exponential distributions. Dierent distributional assumptions, e.g., a dependent joint distribution in which each component is a Weibull random variable might provide a higher "delity model. These distinctions can arise because of important dierences in assumptions made on the underlying physical model governing component failure. Hypothesizing competing models for an underlying phenomenon and understanding the domain of applicability of such models is, of course, central to scienti"c investigation. 5 Conclusions and Some Emerging Tools for Sensitivity Analysis In this report, we have proposed a 10-step sensitivity analysis procedure that we see as practical for large-scale stochastic computer simulation models. And, we have illustrated these ideas on a simple example of a simulation model for system reliability. None of the steps we propose are new. Rather, we rely on the literatures from decision analysis, econometrics, statistics, and simulation to guide what we have proposed. Tornado diagrams provide a simple means for visualizing the in"uence of a signi"cant number of input parameters on an output variable and for understanding to which input variables the output variable is most sensitive. Eective use of tornado diagrams requires the analyst to specify reasonable and commensurate ranges for a collection of input parameters. Both sensitivity plots and spider plots complement tornado diagrams in that they more easily: (i) depict the nonlinear response of an output variable to changes in an input parameter, and (ii) depict the sampling-based error associated with estimating an output measure. We have advocated careful distinction of Monte Carlo sampling-based errors from uncertainties in input parameters. We model the latter type of uncertainty by placing a (possibly joint) probability distribution on the input parameters, 41
and in this setting we seek to understand the resulting probability distribution on the output by propagating the uncertainty through the nonlinear function represented by the simulation model. We discuss how to propagate uncertainty for univariate and multivariate distributions on input parameters. Finally, we describe regression metamodels, and associated experimental design, for understanding sensitivities and interactions between input parameters. A number of important themes repeat throughout our recommendations. These have included the use of common random numbers in order to reduce variability and smooth output analysis. We have also focused on assessing dierences in input parameters that make a dierence in decisions or qualitative characterizations of the system at hand. The relevant ideas we have discussed in this regard include threshold analyses, indierence zones, and establishing dominance relations. Sensitivity analysis simpli"es signi"cantly when using deterministic simulation models. However, our simulations are stochastic and hence proper characterization of both sampling error and un-certainties on input parameters has been a pervasive theme in our presentation. Our discussion is by no means comprehensive for sensitivity analysis of computer simulation models. Alternatives include computing derivatives, which is often termed local sensitivity analysis; see, e.g., Sobol [19]. Speci"cally, with f (x1 , x2 , . . . , xn ) denoting a single output measure from our simulation model, we could estimate the gradient
f f x f (x1 , . . . , xn ) = ,..., . x1 xn Local sensitivity analysis is of particular interest when attempting to optimize f , but optimization over the inputs is not our focus here. More importantly, if we simply report estimates of the partial f derivatives x i
, i = 1, . . . , n, this can mislead with respect to what input parameters are most important because: (i) a per unit change in the temperature of water at a sump pump may not be commensurate with a per gram/fuel-assembly change in debris mass having penetrated the strainer; f
and, (ii) even if we compute x i xi for commensurate values of xi , a linear approximation may be poor over the range of parameters of interest. It is for the reasons just discussed that we advocate the global sensitivity analysis that we have proposed in steps 4-10. Here, the notion of global is speci"ed by the analyst via ranges associated with the input parameters (as we have done in Table 1) rather than rates of change at a single point. These ranges play a central role in sensitivity plots, tornado diagrams, spider plots, and the experimental designs associated with regression metamodels. (Even though we term these global, it is important to recognize that all but the regression metamodel involve changing one parameter at a time.) Such analyses need not associate a probability distribution with the input parameters. However, in our view it is preferable when such probability distributions can be speci"ed because we can make use of them in UQ plots as well as our other sensitivity analyses. In such cases the speci"c endpoints of the ranges become less important. We have recommended using tornado diagrams as an initial tool for assessing the most important input parameters, and using sensitivity plots, UQ plots, spider plots, and metamodels to enable a richer exploration of model sensitivity. It is possible to employ more sophisticated statistical 42
schemes for screening factors using metamodels [25], including sequential bifurcation screening [24]. However, even when such schemes are advocated, it is acknowledged that such approaches have yet to see signi"cant application in practice [9]. Originally developed for interpolation in geostatistical and spatial sampling, Kriging (see, e.g., [5, 17, 18]) has seen widespread successful application in the context of deterministic simulation models as an alternative to regression metamodels. More recently, Kriging metamodels have begun to see application to stochastic simulation models; see, e.g., [1, 7, 2]. While these approaches are a promising alternative to regression metamodels, there are a number of outstanding research issues that remain to be solved [9]. 43
References [1] Ankenman, B., B. L. Nelson, and J. Staum (2008). Stochastic Kriging for simulation metamod-eling. In S. J. Mason, R. R. Hill, L. Moench, and O. Rose (Eds.), Proceedings of the 2008 Winter Simulation Conference, pp. 362-370. [2] Beers, W. V. and J. P. C. Kleijnen (2003). Kriging for interpolation in random simulation. Journal of the Operational Research Society 54, 255-262. [3] Clemen, R. T. and R. Reilly (2001). Making Hard Decisions with Decision Tools. Duxbury, Paci"c Grove, CA. [4] Commission, N. R. (2011). Regulatory Guide 1.174: An Approach for Using Probabilistic Risk Assessment In Risk-Informed Decisions On Plant-Speci"c Changes to the Licensing Basis, Revision 2, Nuclear Regulatory Commission, Washington, DC. [5] Cressie, N. A. C. (1993). Statistics for Spatial Data. Wiley, New York. [6] Eschenbach, T. G. (1992). Spiderplots versus tornado diagrams for sensitivity analysis. Inter-faces 22, 40-46. [7] Gramacy, R. B. and H. K. H. Lee (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association 103, 1119-1130. [8] Kleijnen, J. P. C. (1995). Veri"cation and validation of simulation models. European Journal of Operational Research 82, 145-162. [9] Kleijnen, J. P. C. (2011). Sensitivity analysis of simulation models. In Wiley Encyclopedia of Operations Research and Management Science. Wiley, New York. [10] Kleijnen, J. P. C. and W. V. Groenendaal (1992). Simulation: A Statistical Perspective. Wiley, Chichester. [11] Kynn, M. (2008). The heuristics and biases bias in expert elicitation. Journal of the Royal Statistical Society: Series A 171, 239-264. [12] Letellier, B., T. Sande, and G. Zigler (2013, November). South Texas Project Risk-Informed GSI-191 Evaluation, Volume 3, CASA Grande Analysis. Technical report, STP-RIGSI191-V03, Revision 2. [13] Morton, D. P., J. Tejada, and A. Zolan (2013, September). South Texas Project Risk-Informed GSI-191 Evaluation: Strati"ed Sampling in Monte Carlo Simulation: Motivation, Design, and Sampling Error. Technical report, STP-RIGSI191-ARAI.03, The University of Texas at Austin. 44
[14] Ogden, N., D. P. Morton, and J. Tejada (2013, June). South Texas Project Risk-Informed GSI-191 Evaluation: Filtration as a Function of Debris Mass on the Strainer: Fitting a Parametric Physics-Based Model. Technical report, SSTP-RIGSI191-V03.06, The University of Texas at Austin. [15] OHagan, A., C. E. Buck, A. Daneshkhah, J. Eiser, P. Garthwaite, D. Jenkinson, J. Oakley, and T. Rakow (2006). Uncertain Judgements: Eliciting Expert Probabilities. Wiley, Chichester. [16] Pan, Y.-A., E. Popova, and D. P. Morton (2013, January). South Texas Project Risk-Informed GSI-191 Evaluation, Volume 3, Modeling and Sampling LOCA Frequency and Break Size. Tech-nical report, STP-RIGSI191-V03.02, Revision 4, The University of Texas at Austin. [17] Sacks, J., W. J. Welch, T. Mitchell, and H. Wynn (1989). Design and analysis of computer experiments. Statistical Science 4, 409-435. [18] Santner, T. J., B. J. Williams, and W. I. Notz (2003). The Design and Analysis of Computer Experiments. Springer-Verlag, New York. [19] Sobol, I. M. (2001). Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation 55, 271-280. [20] Tregoning, R., P. Scott, and A. Csontos (2008, April). Estimating Loss-of-Coolant Accident (LOCA) Frequencies Through the Elicitation Process: Main Report (NUREG-1829). NUREG 1829, NRC, Washington, DC. [21] Tversky, A. and D. Kahneman (1981). Rational choice and the framing of decisions. Sci-ence 211, 453-458. [22] Wake"eld, D. and D. Johnson (2013, January). South Texas Project Risk-Informed GSI-191 Evaluation, Volume 2, Probabilistic Risk Analysis: Determination of Change in Core Damage Frequency and Large Early Release Frequency Due to GSI-191 Issues . Technical report, STP-RIGSI191-VO2, Revision 0, ABSG Consulting Inc. [23] WCAP-16793-NP (2011, January). Evaluation of Long-Term Cooling Considering Particulate, Fibrous and Chemical Debris in the Recirculating Fluid, Revision 2. [24] Xu, J., F. Yang, and H. Wan (2007). Controlled sequential bifurcation for software reliability study. In S. G. Henderson, B. Biller, M.-H. Hsieh, J. Shortle, J. D. Tew, and R. R. Barton (Eds.), Proceedings of the 2007 Winter Simulation Conference, pp. 281-288. [25] Yu, H.-F. (2007). Designing a screening experiment with a reciprocal Weibull degradation rate. Computers and Industrial Engineering 52, 175-191. 45
A Appendix: A Sensitivity Analysis for STP GSI-191 In this appendix we apply the proposed framework to analyze the sensitivity of estimates of a risk measure to changes in input parameters of the CASA Grande simulation model, using STP data. The CASA Grande simulation model has the characteristics we describe in the abstract and in Section 1, and the framework we describe in this report was designed with large-scale stochastic simulation models like CASA Grande in mind. In what follows, we brie"y review the sampling scheme within CASA Grande, and we describe the loose coupling between CASA Grande and the PRA model. We describe how we estimate risk in terms of the contribution to core damage frequency (CDF) from GSI-191, in units of events per calendar year (CY), using estimates of the conditional failure probabilities that CASA Grande provides. That is, our risk measure is the change in core damage frequency (CDF) relative to a base CDF due to non-GSI-191 issues. We present the results of 22 scenarios, where each scenario speci"es the values of the input parameters to CASA Grande and where one of these scenarios contains nominal values for the parameters. This presentation includes a tornado diagram representing changes in all parameters we consider, and further includes a sensitivity plot for one of the key inputs. A.1 Step 1: De"ne the Model We refer to Volume 3 [12] for a discussion of the CASA Grande simulation model. One important aspect of this simulation for the purpose of our analysis here is the fact that a strati"ed sampling estimator is employed, in which the strati"cation is on the initiating frequency. The probability distribution governing the initiating frequency is consistent with percentiles from NUREG-1829 [20] as we describe in [16]. This strati"ed sampling estimator can be thought of as an outer loop of replications when running the simulation model, which we refer to as frequency replications. This outer loop facilitates preservation of the probability distribution for initiating frequency in the sense of the uncertainty quanti"cation plots described in Step 7, and the strati"ed estimator further reduces variance versus a nave Monte Carlo estimator. Within each frequency replication, i.i.d. replications are performed in order to estimate con-ditional failure probabilities for each mode of failure (sump and boron "ber limit) and break size (small, medium, and large), conditioned on the pump state as well as the initiating frequency. A strati"cation with 15 cells is used for the strati"ed estimator with respect to the initiating frequency, and importantly, the sampling in distinct cells of the strata is done independently, unless speci"ed otherwise. A sample size within each cell of the strati"cation is selected, as well as the boundaries of each cell of the strata, as indicated in Table 8. See [13] for background on the strati"ed sampling estimator. The right-most column in Table 8 is based on optimization model (10) in [13]. 46
Table 8: Strati"cation of initiating frequency in terms of quantiles of its distribution function F . The probability mass for each cell is indicated, as is the sample size devoted to i.i.d. replications within each cell. Frequency Cell Cell Probability Number of Replication Lower Limit Upper Limit Mass Stat. Replications 1 F 1 (0.000) F 1 (0.045) 0.045 11 2 F 1 (0.045) F 1 (0.115) 0.070 12 3 F 1 (0.115) F 1 (0.195) 0.080 11 4 F 1 (0.195) F 1 (0.260) 0.065 9 5 F 1 (0.260) F 1 (0.295) 0.035 7 6 F 1 (0.295) F 1 (0.365) 0.070 11 7 F 1 (0.365) F 1 (0.435) 0.070 8 8 F 1 (0.435) F 1 (0.510) 0.075 23 9 F 1 (0.510) F 1 (0.620) 0.110 45 10 F 1 (0.620) F 1 (0.685) 0.065 18 11 F 1 (0.685) F 1 (0.720) 0.035 13 12 F 1 (0.720) F 1 (0.830) 0.110 51 13 F 1 (0.830) F 1 (0.955) 0.125 50 14 F 1 (0.955) F 1 (0.990) 0.035 28 15 F 1 (0.990) F 1 (1.000) 0.010 11 A.2 Step 2: Select Outputs of Interest The change in core damage frequency (CDF), when accounting for GSI-191 processes, is selected as the output of interest for this sensitivity study. A.2.1 Core Damage Frequency The method for estimating CDF combines estimates from the CASA Grande simulation model with coecients, as we explain here, from STPs PRA [22]. The CASA Grande simulation model is used to estimate the conditional probabilities of a sump failure and a boron "ber limit failure at various break sizes (small, medium, and large), when we condition on the initiating frequency of a LOCA and the pump state. From the PRA, we can use a base CDF from non-GSI-191 events, as well as the core damage frequencies associated with a sump or boron "ber limit failure, further conditioned on each permutation of pump state and break size. This is formalized as follows. 47
Indices and Sets: i = 1, . . . , F index for cells stratifying frequency replications j = 1, . . . , Mi index for statistical replications k = 1, . . . , N index for set of pump states Events: SL small LOCA ML medium LOCA LL large LOCA P Sk pumps in state k Fi initiating frequency in cell i S sump failure B boron "ber limit CD core damage Parameters: fSL frequency (events/CY) of a Small LOCA fM L frequency (events/CY) of a Medium LOCA fLL frequency (events/CY) of a Large LOCA P (P Sk ) probability mass of P Sk P (Fi ) probability mass of Fi P (SlLOCA, Fi , P Sk ) estimate of probability of S given LOCA = SL, M L, orLL, Fi , P Sk P (BlLOCA, Fi , P Sk ) estimate of probability of B given LOCA = SL, M L, orLL, Fi , P Sk RBASE non-GSI-191 core damage frequency (events/CY) RCD estimate of core damage frequency (events/CY) The three frequencies, fSL , fM L , and fLL , are taken from the right-most column of Table 4-1 from Volume 2 [22]. In the sensitivity analysis computations that follow, we use F = 15 frequency replications, and for the i-th frequency replication cell, we use Mi statistical replications with the values for Mi given in the right-most column of Table 8. As discussed in Volumes 2 and 3 [22, 12], in general we would consider a total of N = 64 pump states, although a reduced number of bounding pump states are used in actual computation. These 64 pump states include the mostly likely Pump State 1, which has all pumps on all three trains available. In the results we present here, we only consider Pump State 1, which has a probability mass of 0.935 when we consider all 64 pump states. In terms of our notation this means we have N = 1 and P (P S1 ) = 1, eectively eliminating the sum over k. If we were to compute P (P Sk ) more generally, we would do so by normalizing the pump state frequencies in the left-hand column of Table 9-1 in Volume 2 [22] or equivalently the same frequencies in Table 2.2.11 in Volume 3 [12]. In our computations we also take RBASE = 0 so that RCD estimates CDF. That said, we develop the formulas that follow for general values of F , Mi , N , and RBASE . 48
Estimating CDF The estimate for CDF is calculated by summing the probability of each pump state-initiating frequency pair, and multiplying each term by the corresponding PRA frequencies coupled with the conditional failure probabilities estimated by CASA Grande. The formula for the estimator is as follows:
F N RCD = RBASE + P (Fi )P (P Sk ) * (4a) i=1 k=1
fSL
- P (SlSL, Fi , P Sk ) + fSL
- P (BlSL, Fi , P Sk ) (4b)
+fM L
- P (SlM L, Fi , P Sk ) + fM L
- P (BlM L, Fi , P Sk ) (4c)
+fLL
- P (SlLL, Fi , P Sk ) + fLL
- P (BlLL, Fi , P Sk ) . (4d)
Each of the six probability estimates, e.g., P (SlLL, Fi , P Sk ), is formed via a sample mean of i.i.d. observations, e.g., P j (SlLL, Fi , P Sk ), j = 1, . . . , Mi , within the CASA Grande simulation. That is, M 1 i j P (SlLL, Fi , P Sk ) = P (SlLL, Fi , P Sk ). (5) Mi j=1 In this way, equation (5), and its "ve analogs (e.g., for P (SlM L, Fi , P Sk )), are substituted for the corresponding terms in equation (4) to form the estimator RCD . Estimating the Variance of the CDF Estimator We must estimate the variance of the estimator RCD in order to quantify its sampling error. There are a total of 6
- F
- N + 1 terms in equation (4) de"ning RCD , and in order to estimate the variance, we must clarify which of these pairs of terms are independent and which are dependent.
First, we assume the terms P (Fi ) and P (P Sk ), as well as the frequencies from the PRA such as RBASE and fLL , are deterministic. Dependency, or the lack thereof, between pairs of estimators like P (SlSL, Fi , P Sk ) and P (BlSL, Fi , P Sk ) depend on how the simulation is performed. The terms across distinct frequency replications and pump states are independent because independent Monte Carlo samples are drawn within each pump state-initiating frequency pair; however, the six terms within each pump state-initiating frequency pair are dependent. Thus, with V (*) and COV (*) denoting the sample variance and sample covariance operators, we have the following equation for V (RCD ), which has 15 sample covariance terms because we have six pump state-initiating frequency pairs: 49
F N V RCD = [P (P Sk )]2 * [P (Fi )]2 * (6a) i=1 k=1
fSL 2
- V P (SlSL, Fi , P Sk ) + fSL 2
- V P (BlSL, Fi , P Sk ) (6b)
L
- V P (SlM L, Fi , P Sk ) + fM L
- V P (BlM L, Fi , P Sk )
2 2
+fM (6c)
2
+fLL
- V P (SlLL, Fi , P Sk ) + fLL 2
- V P (BlLL, Fi , P Sk ) (6d)
+2
- fSL
- fM L
- COV P (SlSL, Fi , P Sk ), P (SlM L, Fi , P Sk ) (6e)
+2
- fSL
- fLL
- COV P (SlSL, Fi , P Sk ), P (SlLL, Fi , P Sk ) (6f)
+2
- fSL
- fSL
- COV P (SlSL, Fi , P Sk ), P (BlSL, Fi , P Sk ) (6g)
+2
- fSL
- fM L
- COV P (SlSL, Fi , P Sk ), P (BlM L, Fi , P Sk ) (6h)
+2
- fSL
- fLL
- COV P (SlSL, Fi , P Sk ), P (BlLL, Fi , P Sk ) (6i)
+2
- fM L
- fLL
- COV P (SlM L, Fi , P Sk ), P (SlLL, Fi , P Sk ) (6j)
+2
- fM L
- fSL
- COV P (SlM L, Fi , P Sk ), P (BlSL, Fi , P Sk ) (6k)
+2
- fM L
- fM L
- COV P (SlM L, Fi , P Sk ), P (BlM L, Fi , P Sk ) (6l)
+2
- fM L
- fLL
- COV P (SlM L, Fi , P Sk ), P (BlLL, Fi , P Sk ) (6m)
+2
- fLL
- fSL
- COV P (SlLL, Fi , P Sk ), P (BlSL, Fi , P Sk ) (6n)
+2
- fLL
- fM L
- COV P (SlLL, Fi , P Sk ), P (BlM L, Fi , P Sk ) (6o)
+2
- fLL
- fLL
- COV P (SlLL, Fi , P Sk ), P (BlLL, Fi , P Sk ) (6p)
+2
- fSL
- fM L
- COV P (BlSL, Fi , P Sk ), P (BlM L, Fi , P Sk ) (6q)
+2
- fSL
- fLL
- COV P (BlSL, Fi , P Sk ), P (BlLL, Fi , P Sk ) (6r)
+2
- fM L
- fLL
- COV P (BlM L, Fi , P Sk ), P (BlLL, Fi , P Sk ) . (6s)
To illustrate computation of the sample variance terms, V (*), and sample covariance terms, COV (*), in equation (6), we give the formulas for the "rst sample variance term from (6b) and the "rst sample covariance term from (6e):
1 1 Mi
2 V P (SlSL, Fi , P Sk ) = P j (SlSL, Fi , P Sk ) P (SlSL, Fi , P Sk ) (7)
Mi Mi 1 j=1 50
and COV (P (SlSL, Fi , P Sk ), P (SlM L, Fi , P Sk ) = (8a) Mi h X i h i 1 1 P j (SlSL, Fi , P Sk ) P (SlSL, Fi , P Sk )
- P j (SlSL, Fi , P Sk ) P (SlM L, Fi , P Sk ) . (8b)
Mi Mi 1 j=1 Further Remarks on Estimating CDF For this sensitivity analysis study, we estimate RCD conditional on being in the pump state in which all three pumps are operating on all three trains; i.e., the most likely pump state called Pump State 1. This eliminates the sum across pump states in the estimators given above and, because we are computing a conditional probability, has the eect of setting the probability mass associated with Pump State 1 equal to 1. In addition we take RBASE = 0 so that RCD estimates CDF rather than CDF. In what follows, we are interested in dierences in RCD under pairs of scenarios, i.e., under two sets of input parameters to CASA Grande. Speci"cally, we have RCD 0 under the nominal settings
of the parameters, RCD under perturbed settings of the parameters, and our interest lies in the estimator RCD = RCD 0 . When computing this dierence, we use common random numbers RCD across the two scenarios to reduce the variance of RCD , which is computed using a straightforward variant of equation (6) in which the sample variance and sample covariance terms have terms like P j (SlSL, Fi , P Sk ) and P (SlSL, Fi , P Sk ) replaced by P j (SlSL, Fi , P Sk ) and P (SlSL, Fi , P Sk ), respectively. Sometimes it is convenient to report and display the ratio of the these frequencies estimates, rather than their dierence; i.e., we report 1 RCD 0
, (9)
RCD and term this the ratio of risk. We again use common random numbers when computing this estimator, and the sample variance of this is estimated via 2 1 RCD 1 RCD 1 RCD 1 V 0
= 0 V (RCD1
)2 0 COV (RCD 1 0
, RCD )+ 0 0
V (RCD ) . (10) RCD [RCD ] 2 RCD RCD A.3 Step 3: Select Inputs of Interest The following input parameters of interest were selected, via multiple discussions of the STP Tech-nical Team, from a larger collection of candidate parameters. The selection was based on: (i) the uncertainty associated with estimates of the values of the parameters, and (ii) the perceived likeli-hood that bias would, based on the teams deliberations, have the largest eect on the estimates of risk, either in the increasing or decreasing direction. 51
- 1. Amount of Latent Fiber in the Pool; 4 levels
- 2. Boron Fiber Limit in the Core; 4 levels
- 3. Debris Transport Fractions Inside the Zone of In"uence; 3 levels
- 4. Chemical Precipitation Temperature; 2 levels
- 5. Total Failure Fractions for Debris Outside the Zone of In"uence (Unquali"ed Coatings); 2 levels
- 6. Chemical Bump-Up Factor; 2 levels
- 7. Fiber Penetration Function; 2 levels
- 8. Size of Zone of In"uence; 2 levels
- 9. Time to Turn O One Spray Pump; 2 levels
- 10. Time to Hot Leg Injection; 2 levels
- 11. Strainer Buckling Limit; 2 levels
- 12. Water Volume in the Pool; 3 levels
- 13. Debris Densities; 2 levels
- 14. Time-Dependent Temperature Pro"les; 2 levels
- 15. Spray Transport Fractions for Debris Outside the Zone of In"uence (Unquali"ed Coatings);
2 levels A short description of the parameters of interest is provided in the following. See also Sections 2 and 5 of Volume 3 [12] for further discussion. Amount of Latent Fiber in Pool There is an amount of existing dust and dirt in the contain-ment, which is based on plant measurement. The latent "ber is assumed to be in the pool at the start of recirculation. This latent debris is therefore available immediately upon start of recirculation, uniformly mixed in the containment pool. During "ll up, this latent debris is also available to penetrate the sump screen. Boron Fiber Limit The boron "ber limit refers to the assumed success criterion, or threshold of concern where boron precipitation would be assumed to occur for cold leg breaks. The "ber limit comes from the testing performed by the vendor that shows no pressure drop occurs with full chemical eects. The assumption is that all "ber that penetrates through the sump screen is deposited uniformly on the core. 52
Debris Transport Fractions in ZOI This refers to the three-zone ZOI debris size distribution. Each dierent type of insulation has a characteristic ZOI which is divided in three sections to take into account the type of damage (debris size distribution) within each zone. Chemical Precipitation Temperature CASA Grande assumes that, once a thin bed of "ber is formed on the strainer, the chemical precipitation bump up factors are applied when the pool temperature reaches the precipitation temperature, de"ned in input. Total Failure Fraction for Debris Outside the ZOI CASA Grande uses a table of total fail-ure fractions that are applied to the transport logic trees. The fraction of each type ("ber, paint and coatings, and so forth) that passes through areas to the pool are used to under-stand what is in the pool as a function of time during recirculation. The total failure fraction multiplies the total inventory of unquali"ed coatings. Chemical Bump Up Factor The chemical bump up factor is used as a multiplier on the con-ventional head loss calculated in CASA Grande. The multiplier is applied if a thin bed is formed and the pool temperature is at or below the precipitation temperature. Fiber Penetration Function The amount of "ber that bypasses the ECCS sump screen (as a fraction) is correlated to the arrival time and the amount of "ber on the screen. The coecients of the correlation de"ne the fractional penetration amounts. Size of ZOI The zone of in"uence (ZOI) is de"ned as a direct function (multiplier) of break size (and nominal pipe diameter). For example, for NUKON "ber, the ZOI (for STP) is 17 times the break diameter. The ZOI is assumed to be spherical unless it is associated with less than a full (double-ended) break in which case it is hemispherical. Otherwise, it is truncated by any concrete walls within the ZOI. Time to Turn O One Spray Pump If three spray pumps start, then by procedure one is se-cured. The time to secure the pump is governed by the operator acting on the conditional action step in the procedure. Time to Hot Leg Injection Similar to the spray pump turn o time, the time to switch one or more trains to hot leg injection operation is governed by procedure. Strainer Buckling Limit The strainer buckling limit is the dierential pressure across the ECCS strainer at which the strainer is assumed to fail mechanically. This limit is based on engi-neering calculations that incorporate a factor of safety. Water Volume in the Pool Depending on the break size, the amount of water that is in the pool, as opposed to held up in the RCS and other areas in containment, is variable. Smaller breaks tend to result in less pool volume than larger breaks. 53
Debris Densities The debris density depends on the amount and type of debris that arrives in the pool. These densities are used in head loss correlations to calculate, for example, debris volume. Time Dependent Temperature Pro"les The temperature of the water in the sump aects air release and vaporization during recirculation. The time-dependent temperature pro"le comes from the coupled RELAP5-3D and MELCOR simulations depending on break size. Spray Transport Fractions for Debris Outside the ZOI CASA Grande uses a table of fail-ure fractions that are applied to the transport logic trees. The fraction of each type ("ber, paint and coatings, and so forth) that passes through areas to the pool are used to understand what is in the pool as a function of time during recirculation. The spray transport fraction is the fraction of failed coatings that wash to the pool during spray operation. The "rst two inputs, latent "ber and boron "ber limit, have four levels, including the nominal case. Items 3 and 12, debris transport within the ZOI and water volume, have three levels. And, all other of the other input parameters have two levels, the nominal case and a perturbation in a single direction. Thus, the number of runs needed to conduct one-at-a-time sensitivity analysis is: 2 * (4 1) + 2 * (3 1) + 11 * (2 1) + 1 = 22. The number of runs needed to conduct a single replicate full factorial design would be 42
- 32
- 211 = 294, 912. We conduct the former, but we do not attempt the latter here.
A.4 Step 4: Choose Nominal Values and Ranges for Inputs The nominal value for an input parameter is sometimes based on the STP Technical Teams best point estimate for that input. However, we sometimes instead select an appropriately conservative nominal value, as is the case for the strainer buckling limit, as we mention above. When selecting ranges for the parameters, sometimes we both increase and decrease a parameter from its nominal value, but other times we only change the parameter in a single direction. We can make the latter choice because we wish to limit changes to directions that we know will increase risk or because the nominal value is already seen as being conservative. Here, we list the 22 scenarios we use for sensitivity analysis. Scenario 0: Nominal-value Case All inputs have a nominal value. Those nominal values are as follows:
- 1. The amount of latent "ber in the pool is 12.5 ft3 .
- 2. The boron "ber limit in the core is 7.5 g/F A.
- 3. The debris transport fractions for debris generated inside the zone of in"uence are given in Table 9.
54
Table 9: Nominal debris transport fractions for debris generated inside the zone of in"uence. Debris Transport Model / LDFG LDFG LDFG Microtherm Qual Coat Crud Debris Type Fines Small Large Fines Fines Fines Blowdown Upper 0.70 0.60 0.22 0.70 0.70 0.70 Blowdown Lower 0.30 0.25 0.00 0.30 0.30 0.30 Washdown Inside 0.53 0.27 0.00 0.53 0.53 0.53 Washdown Annulus 0.47 0.19 0.00 0.47 0.47 0.47 Washdown BC Inside 0.00 0.27 0.00 0.00 0.00 0.00 Washdown BC Annulus 0.00 0.00 0.00 0.00 0.00 0.00 Pool Fill Sump 0.02 0.00 0.00 0.02 0.02 0.02 Pool Fill Inactive 0.05 0.00 0.00 0.05 0.05 0.05 Recirculation Lower 1.00 0.64 0.00 1.00 1.00 1.00 Recirculation Inside 1.00 0.64 0.00 1.00 1.00 1.00 Recirculation Annulus 1.00 0.58 0.00 1.00 1.00 1.00 Erosion Spray 0.00 0.01 0.01 0.00 0.00 0.00 Erosion Pool 0.00 0.07 0.07 0.00 0.00 0.00
- 4. The chemical precipitation temperature is 140 F.
- 5. The total failure fractions for debris generated outside the zone of in"uence (unquali"ed coatings) are given in row 1 of Table 10.
Table 10: Nominal debris failure fractions for debris generated outside the zone of in"uence (un-quali"ed coatings). Debris Transport Model / Epoxy Epoxy Epoxy Epoxy Epoxy Alkyd Baked IOZ Debris Type Fines Fine Small Large Curls Enamel Fines Chips Chips Chips Total Failure Fraction 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 Upper Containment 0.15 0.15 0.15 0.15 0.15 0.54 0.00 0.83 Lower Containment 0.02 0.02 0.02 0.02 0.02 0.46 1.00 0.17 Reactor Cavity 0.83 0.83 0.83 0.83 0.83 0.00 0.00 0.00 Prior to Securing Sprays 0.06 0.06 0.06 0.06 0.06 0.06 0.00 0.06 Recirculation Lower Containment 1.00 0.41 0.00 0.00 1.00 1.00 1.00 1.00 Recirculation Reactor Cavity 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
- 6. The chemical bump-up factors as a function of break size are:
Small Breaks: Truncated Exponential Distribution Mean = 1.25 Minimum = 1 Maximum = 15.3 Medium Breaks: Truncated Exponential Distribution Mean = 1.50 Minimum = 1 Maximum = 18.2 Large Breaks: Truncated Exponential Distribution 55
Mean = 2.00 Minimum = 1 Maximum = 24.0
- 7. The "ber penetration function parameters are:
Fraction of Sheddable Debris: Uniform Distribution Min = 0.00956 Max = 0.02720 Shedding Rate (1/min): Uniform Distribution Min = 0.008236 Max = 0.054600 Eciency Per Gram of Debris: Uniform Distribution Min = 0.000339 Max = 0.003723 Fit Cut Point (g): Uniform Distribution Min = 790 Max = 880 Initial Eciency: Uniform Distribution Min = 0.656 Max = 0.706 Exponential Rate Constant (1/g): Continuous Empirical Distribution Parameters: 0.0011254; 0.10 0.0013078; 0.45 0.0317870; 0.10
- 8. The size of zones of in"uence (R/D) as a function of break size are:
NUKON: 17.0 NUKON 2: 17.0 Microtherm: 28.6 RMI: 1.0 Lead: 1.0 Thermal Wrap: 17.0 IOZ: 1.0 Alkyd: 1.0
- 9. The times to turn o one spray pump (minutes) as a function of break size are:
Small Breaks: 0.0 Medium Breaks: Normal Distribution 56
Mean = 20.0 Standard Deviation = 5.0 Large Breaks: Normal Distribution Mean = 20.0 Standard Deviation = 5.0
- 10. The times to hot leg injection (minutes) as a function of break size are:
Small Breaks: Uniform Distribution Min = 345 Max = 360 Medium Breaks: Uniform Distribution Min = 345 Max = 360 Large Breaks: Uniform Distribution Min = 345 Max = 360
- 11. The strainer buckling limit is 9.35 f t H2 O.
- 12. The water volumes in the pool (ft3 ) as a function of break size are:
Small Breaks: Uniform Distribution Min = 43, 464 Max = 61, 993 Medium Breaks: Uniform Distribution Min = 39, 533 Max = 69, 444 Large Breaks: Uniform Distribution Min = 45, 201 Max = 69, 263
- 13. The debris densities (lbm /f t3 ) are:
LDFG Fines: 2.4 LDFG Small: 2.4 LDFG Large: 2.4 Microtherm Filaments: 2.4
- 14. The temperature pro"les ( F) as a function of break size are given in Figure 22.
57
Temperature Profiles 200 Small & Medium Breaks Large Breaks 175 Temperature (°F) 150 125 100 75 0 100 200 300 400 500 600 700 Time After Break (Hours) (a) Full 700 Hour Pro"le Temperature Profiles 200 Small & Medium Breaks Large Breaks 175 Temperature (°F) 150 125 100 75 0 5 10 15 20 25 30 35 Time After Break (Hours) (b) 36 Hour Pro"le Seen By CASA Grande Figure 22: Temperature pro"les ( F) for small and large breaks. Panel (a) shows the full 700 hour temperature pro"le, while panel (b) shows only the "rst 36 hours, which is the standard run length for CASA Grande scenarios.
- 15. The spray transport fractions for debris generated outside the zone of in"uence (unquali"ed coatings) are given in the "fth row of Table 10.
We use an upper limit of 15 g/F A for the boron "ber limit in Scenario 5 because WCAP-16793 [23] established a boron "ber limit of 15 g/F A as the maximum value for the STP fuel design based on cooling. The reason we decrease below that value in the nominal case, and further in Scenario 4, is to understand the sensitivity to this value since it is not the value for core cooling but rather the value we use for boric acid precipitation failure. At 15 g/F A, there is eectively no resistance 58
to "ow as measured in WCAP-16793. In contrast to the situation for the boron "ber limit, for the strainer buckling limit, the nominal value is conservative, containing a signi"cant safety margin. Hence, we do not make this limit even weaker in our sensitivity analysis, and we only consider an increase in the limit in Scenario 16. Scenario 1: Decrease Latent Fiber The amount of latent "ber in the pool is changed to 6.25 ft3 , whereas the nominal value is 12.5 ft3 . Scenario 2: Increase Latent Fiber I The amount of latent "ber in the pool is changed to 25.0 ft3 , whereas the nominal value is 12.5 ft3 . Scenario 3: Increase Latent Fiber II The amount of latent "ber in the pool is changed to 50.0 ft3 , whereas the nominal value is 12.5 ft3 . Scenario 4: Decrease Boron Fiber Limit The boron "ber limit in the core is changed to 4.0 g/F A, whereas the nominal value is 7.5 g/F A. Scenario 5: Increase Boron Fiber Limit I The boron "ber limit in the core is changed to 15.0 g/F A, whereas the nominal value is 7.5 g/F A. Scenario 6: Increase Boron Fiber Limit II The boron "ber limit in the core is changed to 50.0 g/F A, whereas the nominal value is 7.5 g/F A. Scenario 7: Increase Debris Transport Inside the Zone of In"uence The debris transport fractions for this scenario are given in Table 11. Note that the only increased debris transport fractions are the blowdown and washdown transport fractions for LDFG "nes and small. 59
Table 11: Increased debris transport fractions for debris generated inside the zone of in"uence. Debris Transport Model / LDFG LDFG LDFG Microtherm Qual Coat Crud Debris Type Fines Small Large Fines Fines Fines Blowdown Upper 1.00 1.00 0.22 0.70 0.70 0.70 Blowdown Lower 1.00 1.00 0.00 0.30 0.30 0.30 Washdown Inside 1.00 1.00 0.00 0.53 0.53 0.53 Washdown Annulus 1.00 1.00 0.00 0.47 0.47 0.47 Washdown BC Inside 1.00 1.00 0.00 0.00 0.00 0.00 Washdown BC Annulus 0.00 0.00 0.00 0.00 0.00 0.00 Pool Fill Sump 0.02 0.00 0.00 0.02 0.02 0.02 Pool Fill Inactive 0.05 0.00 0.00 0.05 0.05 0.05 Recirculation Lower 1.00 0.64 0.00 1.00 1.00 1.00 Recirculation Inside 1.00 0.64 0.00 1.00 1.00 1.00 Recirculation Annulus 1.00 0.58 0.00 1.00 1.00 1.00 Erosion Spray 0.00 0.01 0.01 0.00 0.00 0.00 Erosion Pool 0.00 0.07 0.07 0.00 0.00 0.00 Scenario 8: Decrease Debris Transport Inside the Zone of In"uence The debris transport fractions for this scenario are given in Table 12. Note that the only decreased debris transport fractions are the blowdown and washdown transport fraction for LDFG "nes and small. Table 12: Decreased debris transport fractions for debris generated inside the zone of in"uence. Debris Transport Model / LDFG LDFG LDFG Microtherm Qual Coat Crud Debris Type Fines Small Large Fines Fines Fines Blowdown Upper 0.60 0.51 0.22 0.70 0.70 0.70 Blowdown Lower 0.26 0.21 0.00 0.30 0.30 0.30 Washdown Inside 0.45 0.22 0.00 0.53 0.53 0.53 Washdown Annulus 0.40 0.16 0.00 0.47 0.47 0.47 Washdown BC Inside 0.00 0.23 0.00 0.00 0.00 0.00 Washdown BC Annulus 0.00 0.00 0.00 0.00 0.00 0.00 Pool Fill Sump 0.02 0.00 0.00 0.02 0.02 0.02 Pool Fill Inactive 0.05 0.00 0.00 0.05 0.05 0.05 Recirculation Lower 1.00 0.64 0.00 1.00 1.00 1.00 Recirculation Inside 1.00 0.64 0.00 1.00 1.00 1.00 Recirculation Annulus 1.00 0.58 0.00 1.00 1.00 1.00 Erosion Spray 0.00 0.01 0.01 0.00 0.00 0.00 Erosion Pool 0.00 0.07 0.07 0.00 0.00 0.00 Scenario 9: Increase Chemical Precipitation Temperature The chemical precipitation temperature is changed to 160 F, whereas the nominal value is 140 F. 60
Scenario 10: Decrease Total Failure Fractions for Debris Outside the Zone of In"uence (Unquali"ed Coatings) Table 13: Decreased total failure fractions for debris generated outside the zone of in"uence (un-quali"ed coatings). Debris Transport Model / Epoxy Epoxy Epoxy Epoxy Epoxy Alkyd Baked IOZ Debris Type Fines Fine Small Large Curls Enamel Fines Chips Chips Chips Total Failure Fraction 0.80 0.80 0.80 0.80 0.80 0.43 0.43 0.92 Upper Containment 0.15 0.15 0.15 0.15 0.15 0.54 0.00 0.83 Lower Containment 0.02 0.02 0.02 0.02 0.02 0.46 1.00 0.17 Reactor Cavity 0.83 0.83 0.83 0.83 0.83 0.00 0.00 0.00 Prior to Securing Sprays 0.06 0.06 0.06 0.06 0.06 0.06 0.00 0.06 Recirculation Lower Containment 1.00 0.41 0.00 0.00 1.00 1.00 1.00 1.00 Recirculation Reactor Cavity 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Scenario 11: Increase Chemical Bump-Up Factor The chemical bump-up factors as a function of break size are changed as follows: Small Breaks: Truncated Exponential Distribution Mean = 1.875 Minimum = 1 Maximum = 30.0 Medium Breaks: Truncated Exponential Distribution Mean = 2.25 Minimum = 1 Maximum = 30.0 Large Breaks: Truncated Exponential Distribution Mean = 3.00 Minimum = 1 Maximum = 30.0 Scenario 12: Increase Fiber Penetration (Lower Envelope) The "ber penetration function parameters are no longer sampled from distributions. Instead, they are now constants with the following values: Fraction of Sheddable Debris: 0.0196 Shedding Rate (1/min): 0.0538 Eciency Per Gram of Debris: 0.0003391 Fit Cut Point (g): 880 Initial Eciency: 0.656 61
Exponential Rate Constant (1/g): 0.0013 Scenario 13: Decrease Size of Zone of In"uence (ZOI) The size of zones of in"uence (R/D) as a function of break size are changed as follows: NUKON: 12.75 NUKON 2: 12.75 Microtherm: 21.45 RMI: 0.75 Lead: 0.75 Thermal Wrap: 12.75 IOZ: 0.75 Alkyd: 0.75 Scenario 14: Increase Time to Turn O One Spray Pump The times to turn o one spray pump (minutes) as a function of break size are changed as follows: Small Breaks: 0.0 Medium Breaks: Normal Distribution Mean = 1440.0 Standard Deviation = 5.0 Large Breaks: Normal Distribution Mean = 1440.0 Standard Deviation = 5.0 Scenario 15: Increase Time to Hot Leg Injection The times to hot leg injection (minutes) as a function of break size are changed as follows: Small Breaks: 450 Medium Breaks: 450 Large Breaks: 450 Scenario 16: Increase Strainer Buckling Limit The strainer buckling limit is changed to 10.30 f t H2 O, whereas the nominal value is 9.35 f t H2 O. 62
Scenario 17: Decrease Water Volume The water volumes in the pool (ft3 ) as a function of break size are changed as follows: Small Breaks: Uniform Distribution Min = 39, 191 Max = 56, 720 Medium Breaks: Uniform Distribution Min = 34, 084 Max = 63, 995 Large Breaks: Uniform Distribution Min = 39, 478 Max = 63, 540 Scenario 18: Increase Water Volume The water volumes in the pool (ft3 ) as a function of break size are changed as follows: Small Breaks: Uniform Distribution Min = 48, 737 Max = 67, 266 Medium Breaks: Uniform Distribution Min = 44, 982 Max = 74, 893 Large Breaks: Uniform Distribution Min = 50, 924 Max = 74, 986 Scenario 19: Increase Debris Densities The debris densities (lbm /f t3 ) are: LDFG Fines: 3.0 LDFG Small: 3.0 LDFG Large: 3.0 Microtherm Filaments: 3.0 Scenario 20: Decrease Time-Dependent Temperature Pro"les The modi"ed temperature pro"les ( F) as a function of break size are given in Figure 23. 63
Temperature Profiles 200 Small & Medium Breaks Large Breaks 175 Temperature (°F) Small & Medium Breaks (Low) Large Breaks (Low) 150 125 100 75 0 100 200 300 400 500 600 700 Time After Break (Hours) (a) Full 700 Hour Pro"le Temperature Profiles 200 Small & Medium Breaks Large Breaks 175 Temperature (°F) Small & Medium Breaks (Low) Large Breaks (Low) 150 125 100 75 0 5 10 15 20 25 30 35 Time After Break (Hours) (b) 36 Hour Pro"le Seen By CASA Grande Figure 23: Temperature pro"les ( F) for small and large breaks. Panel (a) shows the full 700 hour temperature pro"le, while panel (b) shows only the "rst 36 hours, which is the standard run length for CASA Grande scenarios. 64
Scenario 21: Increase Spray Transport Fraction for Debris Outside the Zone of In"u-ence (Unquali"ed Coatings) Table 14: Increased spray transport fractions for debris generated outside the zone of in"uence (unquali"ed coatings). Debris Transport Model / Epoxy Epoxy Epoxy Epoxy Epoxy Alkyd Baked IOZ Debris Type Fines Fine Small Large Curls Enamel Fines Chips Chips Chips Total Failure Fraction 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 Upper Containment 0.15 0.15 0.15 0.15 0.15 0.54 0.00 0.83 Lower Containment 0.02 0.02 0.02 0.02 0.02 0.46 1.00 0.17 Reactor Cavity 0.83 0.83 0.83 0.83 0.83 0.00 0.00 0.00 Prior to Securing Sprays 0.12 0.12 0.12 0.12 0.12 0.12 0.00 0.12 Recirculation Lower Containment 1.00 0.41 0.00 0.00 1.00 1.00 1.00 1.00 Recirculation Reactor Cavity 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 A.5 Step 5: Estimating Model Outputs under Nominal Values of Input Parameters Table 15 presents the results from running CASA Grande with all parameters set at their nominal values (Scenario 0), using the formula for RCD from equation (4), along with estimates of sampling error. We see the point estimate of risk in terms of CDF (events/CY) is 1.817E-08, with lower and upper 95% con"dence limits of 1.626E-08 and 2.009E-08, respectively. If we take 1.00E-06 as a threshold of interest then our estimate of RCD is not close to this threshold when all input parameters are set to their nominal values. Table 15: Results for Scenario 0 with all input parameters set to their nominal values. The "rst column is RCD in units of events/CY while the remaining columns characterize the sampling error. Mean 95% CI 95% CI 95% CI CI HW % Risk Half-Width Lower Limit Upper Limit of Mean 1.817E-08 1.914E-09 1.626E-08 2.009E-08 10.53% A.6 Step 6: One-Way Sensitivity Analysis: Sensitivity Plots and Tornado Diagrams In this step, we present results associated with running all 22 scenarios both in tabular form and using a tornado diagram, and we further present a one-way sensitivity plot. We begin with numerical results in tabular form for both the absolute risk (again, conditional on being in the state with all pumps working in all three trains) and for the dierence in CDF with respect to a perturbation of the input parameters relative to the nominal parameter values. Then, we present a tornado diagram corresponding to changing the 15 inputs as we detail in Section A.4. Finally, we present a one-way sensitivity plot for the boron "ber limit, which, after examining the tornado plot and tables of results, appears to be the input parameter to which our estimate of CDF is most sensitive. 65
Tables of Results Table 16 presents the numerical results for each of the 22 scenarios. In this table, we present our estimate of the CDF (Mean Risk) associated with each scenario, and we provide 95% con"dence interval (CI) limits associated with this estimate. We also present the ratio of the 95% CI half-width to the mean as a percentage. Table 16: CDF estimates, and sampling error, for all scenarios. The "rst two columns provide the scenario number and the parameter being changed; see Section A.4. The third column reports how we anticipate the CDF will change, and the remaining columns read in the same manner as in Table 15. Sensitivity Expected Mean 95% CI 95% CI 95% CI CI HW % Measure Direction Risk Half-Width LL UL of Mean 0 Baseline None 1.817E-08 1.914E-09 1.626E-08 2.009E-08 10.53% 1 Latent Fiber Low (6.25 ft3 ) Decrease 1.905E-08 1.928E-09 1.712E-08 2.098E-08 10.12% 2 Latent Fiber High (25 ft3 ) Increase 1.669E-08 1.770E-09 1.492E-08 1.846E-08 10.61% 3 Latent Fiber Very High (50 ft3 ) Increase 3.394E-08 1.447E-08 1.947E-08 4.840E-08 42.63% 4 Boron Fuel Limit (4.0 g/F A) Increase 1.690E-06 1.146E-06 5.445E-07 2.836E-06 67.79% 5 Boron Fuel Limit (50 g/F A) Decrease 1.308E-08 1.412E-09 1.167E-08 1.449E-08 10.80% 6 Boron Fuel Limit (15 g/F A) Decrease 1.329E-08 1.415E-09 1.188E-08 1.471E-08 10.65% 7 Debris Transport Inside ZOI High Increase 7.896E-08 2.250E-08 5.645E-08 1.015E-07 28.50% 8 Debris Transport Inside ZOI Low Decrease 1.241E-08 1.493E-09 1.092E-08 1.390E-08 12.03% 9 Chemical Temp High Increase 1.905E-08 1.937E-09 1.712E-08 2.099E-08 10.17% 10 Total Failure % Outside ZOI Low (80%) Decrease 1.770E-08 1.878E-09 1.582E-08 1.958E-08 10.61% 11 Bump Factor High Increase 2.287E-08 2.024E-09 2.085E-08 2.490E-08 8.85% 12 Penetration Low Envelope Increase 1.552E-07 1.696E-08 1.382E-07 1.721E-07 10.93% 13 ZOI Size Small Decrease 6.795E-09 8.275E-10 5.967E-09 7.622E-09 12.18% 14 Turn O 1 Spray Longer Decrease 1.569E-08 1.763E-09 1.393E-08 1.745E-08 11.23% 15 Hot Leg Injection Longer Increase 1.962E-08 1.954E-09 1.766E-08 2.157E-08 9.96% 16 Strainer Limit Higher Decrease 1.639E-08 1.801E-09 1.459E-08 1.819E-08 10.99% 17 Water Volume Low Increase 2.001E-08 2.027E-09 1.798E-08 2.203E-08 10.13% 18 Water Volume High Decrease 1.655E-08 1.776E-09 1.477E-08 1.833E-08 10.73% 19 Debris Density High Increase 2.567E-08 2.353E-09 2.331E-08 2.802E-08 9.17% 20 Temperature Pro"les Low Increase 1.963E-08 1.991E-09 1.764E-08 2.162E-08 10.14% 21 Spray Transport % Outside ZOI High (12%) Increase 1.798E-08 1.914E-09 1.606E-08 1.989E-08 10.65% Table 17 focuses on the dierences and the ratios of the CDF under the nominal parameter values and under the perturbed parameter values; i.e., we report RCD /R0 CD in the third column
(Ratio) and RCD = RCD RCD in the fourth column (Mean Di), where RCD 0 is the point
estimate of CDF under the nominal scenario and RCD is that under the perturbation scenarios. 66
Table 17: CDF estimates for all scenarios in comparison with the nominal case. This comparison is performed via the ratio and the dierence. The CI statements are for the dierence, and the "nal column indicates whether the dierence is statistically signi"cant at a 95% con"dence level. Sensitivity Ratio Mean 95% CI 95% CI 95% CI Sig Measure Di Half-Width LL UL Di ? 0 Baseline 1.00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 No 1 Latent Fiber Low (6.25 ft3 ) 1.05 8.76E-10 5.37E-10 3.39E-10 1.41E-09 Yes 2 Latent Fiber High (25 ft3 ) 0.92 -1.49E-09 8.00E-10 -2.29E-09 -6.87E-10 Yes 3 Latent Fiber Very High (50 ft3 ) 1.87 1.58E-08 1.45E-08 1.27E-09 3.02E-08 Yes 4 Boron Fuel Limit (4.0 g/F A) 93.01 1.67E-06 1.15E-06 5.24E-07 2.82E-06 Yes 5 Boron Fuel Limit (50 g/F A) 0.72 -5.10E-09 1.09E-09 -6.18E-09 -4.01E-09 Yes 6 Boron Fuel Limit (15 g/F A) 0.73 -4.88E-09 1.08E-09 -5.96E-09 -3.80E-09 Yes 7 Debris Transport Inside ZOI High 4.34 6.08E-08 2.22E-08 3.86E-08 8.30E-08 Yes 8 Debris Transport Inside ZOI Low 0.68 -5.76E-09 1.05E-09 -6.81E-09 -4.72E-09 Yes 9 Chemical Temp High 1.05 8.78E-10 4.22E-10 4.56E-10 1.30E-09 Yes 10 Total Failure % Outside ZOI Low (80%) 0.97 -4.73E-10 3.07E-10 -7.80E-10 -1.66E-10 Yes 11 Bump Factor High 1.26 4.70E-09 9.52E-10 3.75E-09 5.65E-09 Yes 12 Penetration Low Envelope 8.54 1.37E-07 1.69E-08 1.20E-07 1.54E-07 Yes 13 ZOI Size Small 0.37 -1.14E-08 1.52E-09 -1.29E-08 -9.86E-09 Yes 14 Turn O 1 Spray Longer 0.86 -2.49E-09 7.02E-10 -3.19E-09 -1.78E-09 Yes 15 Hot Leg Injection Longer 1.08 1.44E-09 2.46E-09 -1.02E-09 3.90E-09 No 16 Strainer Limit Higher 0.90 -1.78E-09 6.57E-10 -2.44E-09 -1.13E-09 Yes 17 Water Volume Low 1.10 1.83E-09 5.45E-10 1.29E-09 2.38E-09 Yes 18 Water Volume High 0.91 -1.62E-09 6.87E-10 -2.31E-09 -9.38E-10 Yes 19 Debris Density High 1.41 7.49E-09 1.36E-09 6.14E-09 8.85E-09 Yes 20 Temperature Pro"les Low 1.08 1.46E-09 5.51E-10 9.05E-10 2.01E-09 Yes 21 Spray Transport % Outside ZOI High (12%) 0.99 -1.96E-10 2.08E-10 -4.04E-10 1.25E-11 No Tornado Diagram and Analysis Because the CD frequencies we estimate are so small, it is useful to present a tornado diagram and one-way sensitivity plot for the ratios on a logarithmic scale; see also our discussion surrounding Table 17. Figure 24 is a tornado diagram for the 15 input parameters we varied for this sensitivity study. Because we report results for the ratio RCD /R0 , if the ratio has value 10, it means that CD the point estimate of the CDF under the perturbed scenario is 10 times greater than that under the nominal scenario. The CI bounds for the risk ratios are calculated using equation (10). In what follows we examine six factors to which the estimate of CDF seems to have the greatest sensitivity. Changes in the other input parameters lead to more modest changes in the CDF estimate, with percentage dierences of less than 30%. 67
Tornado Diagram: Risk Ratio of Risk under the Scenarios to Risk under Nominal Parameter Values Decreasing Risk Increasing Risk 0.10 1.00 10.00 100.00 1,000.00 Boron Fuel Limit (4.0 g/FA - 50 g/FA) Penetration Low Envelope Debris Transport Inside ZOI Decreased Parameter Values ZOI Size Small Increased Parameter Values Scenario Descriptions Latent Fiber (6.25 ft^3 - 50 ft^3) Debris Density High Bump Factor High Water Volume Turn Off 1 Spray Longer Strainer Limit Higher Temperature Profiles Low Hot Leg Injection Longer Chemical Temp High Total Failure % Outside ZOI Low (80%) Spray Transport % Outside ZOI High (12%) Figure 24: Tornado diagram in log space for ratios of the CDF estimates (risk); i.e., we plot
/R0 , along with the corresponding con"dence intervals for each endpoint of the horizontal RCD CD bars. Values of the ratio that exceed one correspond to an increase in risk relative to the nom-inal case. The scenario numbers in the "rst column of Tables 16 and 17 map the brief scenario descriptions here to the richer scenario descriptions in Section A.4.
We see from the tornado diagram in Figure 24 that the boron "ber limit appears to be the factor to which our estimate of CDF is most sensitive. Increasing the limit from its nominal value of 7.5 g/F A should decrease the CDF, and decreasing the limit should have the opposite eect. This holds because under a larger limit, more "ber can penetrate the strainer without the simulation model declaring a failure. Increasing the limit from 7.5 g/F A to 15.0 g/F A decreases CDF by about 27%. Increasing the value further to 50.0 g/F A leads to little further decrease in CDF. However, decreasing the value from 7.5 g/F A to 4.0 g/F A increases the point estimate of CDF to 1.69E-06, larger than the nominal point estimate by a factor of 93. The tornado diagram suggests the next perturbation to which the estimate of CDF is most sensitive involves changing the "ltration function to its lower envelope using the estimates provided in [14] and reported in Section A.4. The "ltration function aects how much "ber penetrates the strainer. Modifying the function so that less mass is "ltered means that more mass penetrates, and hence, we anticipate this will increase CDF. At this lower envelope, the CDF estimate increases by a factor of 8.5 to 1.55E-07. The transport matrix for debris inside the ZOI governs the amount of each type of debris 68
transported to the sump. When this matrix has smaller transport fractions, we would expect the estimated CDF to decrease because less debris reaches the strainer, and larger transport fractions should have the opposite eect. Under the perturbation to smaller transport fractions speci"ed in Section A.4, the point estimate of CDF decreases by about 32% to 1.24E-08, and with more debris transported to the strainer the CDF estimate grows by a factor of 4.34 to 7.9E-08. Next, we examine the eect of the size of the ZOI. We regard our nominal estimates of the ZOI size as conservative (larger ZOI than expected), and so we examine the result of reducing the size of the ZOI. A smaller ZOI means less debris will be generated, which should reduce CDF. Reducing the size of the ZOI as speci"ed in Section A.4, decreases CDF by about 63% to 6.8E-09. We now examine the eect of latent "ber in the sump. We anticipate that increasing the amount of latent "ber in the sump will increase CDF as more "ber reaches the strainer and can penetrate to the core. As the tornado diagram and Table 17 indicate, increasing the amount of latent "ber from 12.5 ft3 to 50.0 ft3 leads to an 87% increase in CDF. However, more modest changes in latent "ber produce counterintuitive results, as shown in Table 17. A decrease in the amount of latent "ber from 12.5 ft3 to 6.25 ft3 leads to a 5% increase in the CDF estimate, and an increase in the amount of latent "ber from 12.5 ft3 to 25.0 ft3 leads to a 8% reduction in the CDF estimate. The reason for these counterintuitive results is as follows. CASA Grande assumes that some fraction of latent "ber is deposited on the screen when the simulation model is initialized. This latent "ber is not eligible to penetrate the screen, although it is eligible to penetrate by shedding. Of course, in reality some of this "ber will penetrate the screen, and hence the STP Technical Team may suggest a minor modi"cation to CASA Grande in this vein. Increasing the debris densities for low-density "berglass and Microtherm "laments leads to a 41% increase in the CDF estimate to a value of 2.6E-08. Even though the change is small in magnitude, we close this section by discussing a counter-intuitive result. Increasing the spray transport fraction for debris outside the ZOI (unquali"ed coatings failure) counterintuitively results in eectively no change or a slight reduction in debris bed head loss and consequently very little change to CDF. This result is caused by a reduction in the overall (composite) surface-area-to-volume ratio (Sv ) of the debris bed when spray transport applied to failed coatings is increased as a single parameter in Scenario 21. The eect is caused by competing ratios of constituents in the weighted average Sv . Speci"cally, enamel coatings have a large inventory and a small diameter, so the constituent Sv for enamel is large, but enamel coatings were assigned a spray transport fraction of zero. As quantities of other unquali"ed coating types increase with increasing spray transport fraction, the relative proportion of enamel decreases, re-sulting in a lower aggregate Sv . The eect is observed regardless of the weighting scheme chosen for Sv , and both the magnitude and direction of the eect depends on the relative inventories and particle diameters that are speci"ed for the coatings debris types. In treatment of head loss through porous media, the parameter Sv represents the total surface area inside of the bed that can induce drag on the internal "ow. Addition of any material inside of the same bed thickness should both increase Sv and decrease porosity. The unusual dependence of the composite Sv on relative debris 69
quantities and characteristics suggests that traditional use of a particulate-weighted bed property is not appropriate. If a composite parameter is needed, total available drag area should be averaged over the spatial dimensions of the bed and not over aggregate properties of the debris elements themselves. One-Way Sensitivity Plot and Analysis Given the signi"cant eect of the boron "ber limit on the CDF estimate, we explore this sensitivity further via a one-way sensitivity plot. Table 18 contains the results of running the simulation with boron "ber limit values of 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.5, and 8.5 (g/F A). Figure 25 is a one-way sensitivity plot of the CDF estimate versus boron "ber limit over this range. As with the tornado diagram, we present the ratio of the CDF estimate for each of these values to the CDF estimate at the nominal value of 7.5 g/F A, and we present these ratios on a log scale. The CI bounds for the risk ratios are calculated using equation (10). We employ common random numbers in the simulation runs across these dierent "ber limit values. Table 18: CDF (Mean Risk) as a function of the boron "ber limit. Boron Fiber Limit (g/F A) Mean Risk 95% CI HW Ratio 4.0 1.690E-06 1.146E-06 93.01 4.5 5.860E-07 8.359E-07 32.24 5.0 1.059E-07 5.789E-08 5.83 5.5 6.242E-08 3.961E-08 3.43 6.0 3.699E-08 2.038E-08 2.04 6.5 2.050E-08 1.999E-09 1.13 7.0 1.931E-08 1.950E-09 1.06 7.5 1.817E-08 1.914E-09 1.00 8.0 1.700E-08 1.790E-09 0.94 8.5 1.658E-08 1.729E-09 0.91 70
Sensitivity Plot: Common Random Numbers 1000.0 Mean Risk Risk = 1.0E-6 100.0 Risk Ratios Increasing Risk 10.0 1.0 Decreasing Risk 0.1 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 Boron Fuel Limit (g/FA) Figure 25: Sensitivity plot for boron "ber limit (g/F A). From Figure 25 we see there is little change in the CDF estimate as the "ber limit ranges from 6.5 g/F A to 8.5 g/F A. However, the CDF estimate grows quickly as we decrease the "ber limit from 6.5 g/F A. Analysis of Amount of Fiber Penetration Using Perturbed Filtration Function In addition to examining the CDF estimate associated with each of the scenarios described above, another performance measure of interest is the amount of "ber penetrating the core (g/F A) for dierent scenarios. In particular, it is interesting to compare the average amount of "ber penetrating the core using the nominal "ltration function (scenario 0), and using the lower envelope of the "ltration function (scenario 12). In Table 19, we present the results of a statistical analysis on "ber penetration for these two cases. When using the nominal "ltration function, the average amount of "ber penetrating the strainer in CASA Grande is 0.318 g/F A, and when using the lower envelope of the "ltration function, this same value is 0.648 g/F A. We see that approximately twice as much "ber penetrates when using the lower envelope "ltration function. The far right column of Table 19 presents statistical information about the dierence between these two scenarios in terms of "ber penetration. We can see the a 95% con"dence interval on the mean dierence does not include zero, which indicates the dierence is statistically signi"cant, and of course, twice as much "ber penetrating the strainer is also signi"cant from a practical perspective. 71
Table 19: Statistical analysis of "ber penetration under the nominal settings of the parameters and when we use the lower envelope of the "ltration function. Measure/Case Nominal Lower Envelope Dierence Mean 0.318 0.648 0.330 Variance 0.063 0.231 0.289 Standard Deviation 0.251 0.481 0.538 Number of Observations 308 308 308 Con"dence Level 0.95 0.95 0.95 CI Half-Width 0.028 0.054 0.060 CI Lower Limit 0.290 0.595 0.270 CI Upper Limit 0.346 0.702 0.390 p-value - - 3.673E-23 Signi"cant Dierence? - - Yes Summary This appendix has focused on identifying the input parameters to which the performance measure of CDF (change in core damage frequency due to GSI-191 issues) is most sensitive. In general, we estimate CDF by coupling conditional failure probabilities, as estimated by the CASA Grande simulation model, with: (i) the frequency of small, medium, and large LOCA events, and (ii) the probability mass function governing the plant having access to a set of ECCS pumps. In the analysis we presented here, we have assumed that the most likely pump-state case, in which the plant has access to all pumps, occurs with probability one. For this sensitivity analysis, the STP Technical Team selected a total of 15 input parameters to the CASA Grande simulation model. (We note that some parameters actually correspond to a collection of parameters; e.g., we simultaneously change a set of debris transport fractions.) Nominal values for these parameters correspond to the analysis performed in Volumes 2 and 3 [22, 12]. Along with this nominal scenario, 21 further scenarios corresponded to changing the values of these 15 parameters, one at a time. Some parameters were changed in only one direction, and other parameters were both increased and decreased. With the ranges for these parameters in hand, we constructed a tornado diagram characterizing the sensitivity of CDF to changes in the input parameters. Key to our analysis is that the perturbations to these 15 input parameters are commensurate, meaning that they represent changes to the nominal case that have comparable likelihood, as judged by the STP Technical Team. Our CDF estimate is most sensitive to three parameters that concern: (i) how much debris is required to trigger an in-vessel failure (boron "ber limit), (ii) the fraction of debris that penetrates the sump strainer ("ber penetration function), and (iii) the fraction of debris of dierent types that is transported from dierent locations during dierent operational phases (debris transport fractions in ZOI). The eect of the boron "ber limit exceeds that of the next most sensitive parameter by an order of magnitude, and so we examined CDF versus the boron "ber limit in further detail via a 72
one-way sensitivity plot. The growth in CDF is modest as we decrease the boron "ber limit from 7.5 grams per fuel assembly (g/F A) to 6.5 g/F A, but then we see sharp growth in CDF as we further decrease this limit. The appendix of this report applies the initial steps of the sensitivity analysis procedure we propose. Additional analysis will be carried out. We will seek to understand, conditional upon a sump or boron "ber limit failure occurring, which weld locations are most likely to have experienced a break. We will form a spider plot (step 8) for the three or four most sensitive parameters. We will construct a meta-model of the type indicated in step 10 of our framework. And, we intend to include further pump states, beyond the most likely state considered here. 73
A Practical Guide to Sensitivity Analysis of a Large-Scale Computer Simulation Model David Morton Jeremy Tejada Alex Zolan The University of Texas at Austin February 2014
Contents
- Background on applying sensitivity analysis to CASA Grande simulation model
- Practical step-by-step guide to sensitivity analysis: 10-step process
- Illustrative example
- Results from CASA Grande
Background
- Distinguish three sources of error:
- sampling-based errors due to Monte Carlo simulation - errors due to uncertainty in model input - errors due to a lack of fidelity of the model
- Include sampling errors via error bars
- Use UQ plots to characterize second uncertainty
- We do not address model uncertainty here
- Use common random numbers (CRNs) to reduce the variance when comparing performance measures
10-Step Sensitivity Analysis Process
- Step 1: Define the Model
- Step 2: Select Outputs of Interest
- Step 3: Select Inputs of Interest
- Step 4: Choose Nominal Values and Ranges for Inputs
- Step 5: Estimate Model Outputs under Nominal Input Values
- Step 6: One-Way Sensitivity Analysis: Sensitivity Plots & Tornado Diagrams
- Step 7: One-Way Sensitivity Analysis: UQ Plots
- Step 8: One-Way Sensitivity Analysis: Spider Plots
- Step 9: Two-way Sensitivity Analysis: Two-way Sensitivity Plots
- Step 10: Metamodels & Design of Experiments
Step 1: Define the Model
- 1,, 2,, 3, 4 : failure rates
- t0 : desired lifetime of system
- t : time required to perform PM on component 3
- k: PM benefit factor 3/k
Step 2: Select Outputs of Interest
- 9:3/&)(0()#" *#
- ;N?<30()#'" " )0
Step 3: Select Inputs of Interest
- BASE Option:
- 1,, 2,, 3, 4 : failure rates - t0 : desired system lifetime
- PM Option:
- Base option parameters plus: - t : time required to perform PM on component 3 - k : PM benefit factor
- Component 3s failure rate drops to 3 /k, where @
Step 4: Choose Nominal Values and Ranges for Inputs &#! #! & (! & ! &
@6@7#%$) (8 A?? @D? AD?
A6@7#%$) (8 A?? @D? AD?
B6@7#%$) (8 D? AD FD
C6@7#%$) (8 D? AD FD
?7#%$) (8 @G @A AC
7#%$) (8 @ ?4D B
A @ D
Step 5: Estimate Outputs under Nominal Values of Input Parameters
&# &# "&! $ $
;N?< ?4FDAM?4?AE ?4FGDM?4?AD
9: CC4?@MA4B? DC4HFMA4HH
! " &# &# "&!" #
;N?< ?4?BBM?4?@C ?4?BBM?4?BE
9: @?4HDM@4BH @B4@?MB4GA
Step 6: 1-Way Sensitivity Analysis: Tornado Diagrams
Step 6: 1-Way Sensitivity Plots Step 6: 1-Way Sensitivity Plots Step 6: 1-Way Sensitivity Analysis: Tornado Diagrams
Step 7: 1-Way Sensitivity Analysis: Uncertainty Quantification Plots
Step 8: 1-Way Sensitivity Analysis: Spider Plots
Step 8: 1-Way Sensitivity Analysis: Spider Plots
Step 9: Two-way Sensitivity Analysis: Two-way Sensitivity Plots
Step 10: Metamodels & Design of Experiments
- So far, change one or two input parameters at a time
- Metamodel (aka, response surface or surrogate model) built on an experimental design, better captures interaction effects
- A parsimonious metamodel is a polynomial regression model of low degree:
Summary
- Proposed a 10-step sensitivity analysis procedure and illustrated ideas on a simple example
- Recommended using tornado diagrams as initial tool for assessing the input parameters to which output is most sensitive
- Recommended using sensitivity plots, UQ plots, spider plots, and metamodels for a richer exploration of model sensitivity
Appendix Sensitivity Analysis for STP GSI-191
Step 1: Define the Model
- We wont detail CASA Grande here (Volume 3)
- Use CASA Grande to estimate probability of sump failure and boron fiber limit failure, conditional on small, medium & large breaks
- Estimate change in core damage frequency (CDF) in events/year due to GSI-191 issues using these failure probabilities and link to PRA
- All results are conditional on all pumps working
Step 2: Select Outputs of Interest
- Change in core damage frequency (CDF)
- Sometimes, we report ratio of CDF estimate for a scenario to CDF estimate for baseline and call this the risk ratio
- Use stratified sampling on initiating frequency
- Use IID replications within each cell of stratification
- Use common random numbers across scenarios; i.e., use CRNs across specified changes in input parameters
Step 2: Outputs: Estimating CDF Indices and Sets: i = 1, . . . , F index for cells stratifying frequency replications k = 1, . . . , N index for set of pump states Events: SL, M L, LL small, medium, large LOCA P Sk pumps in state k Fi initiating frequency in cell i S sump failure B boron "ber limit failure CD core damage Parameters: fSL , fM L , fLL frequency (events/CY) of a small, medium, large LOCA P (P Sk ) probability mass of P Sk P (Fi ) probability mass of Fi P (SlLOCA, Fi , P Sk ) estimate of probability of S given LOCA = SL, M L, orLL, Fi , P Sk P (BlLOCA, Fi , P Sk ) estimate of probability of B given LOCA = SL, M L, orLL, Fi , P Sk RBASE non-GSI-191 core damage frequency (events/CY) RCD estimate of core damage frequency (events/CY)
Step 2: Outputs: Estimating CDF CDF = RCD RBASE
F N
= P (Fi )P (P Sk )
- i=1 k=1
fSL
- P (SlSL, Fi , P Sk ) + fSL
- P (BlSL, Fi , P Sk )
+fM L
- P (SlM L, Fi , P Sk ) + fM L
- P (BlM L, Fi , P Sk )
+fLL
- P (SlLL, Fi , P Sk ) + fLL
- P (BlLL, Fi , P Sk )
- We report results with:
- fSL , fML , fLL from Volume 2s Table 4-1 - P(all pumps working)=1 - P(Fi ): Bounded Johnson fit to NUREG-1829
- We form a variance estimate for the above estimator
Step 3: Select Inputs of Interest
- Amount of Latent Fiber in Pool: Existing dust/dirt in containment, based on plant measurement. Assumed to be in the pool at start of recirculation, uniformly mixed. During fill up, latent debris available to penetrate sump screen.
- Boron Fiber Limit: Refers to threshold where boron precipitation occurs for cold leg breaks. Fiber limit comes from vendor testing that shows no pressure drop occurs with full chemical effects. Assume all fiber that penetrates sump screen deposits uniformly on core.
- Debris Transport Fractions in ZOI: Refers to three-zone ZOI debris size distribution. Each insulation type has characteristic ZOI divided in three sections to account for type of damage within each zone.
Step 3: Select Inputs of Interest
- Chemical Precipitation Temperature: CASA Grande assumes that, once a thin bed of fiber forms on strainer, chemical precipitation bump up factors apply when pool temperature reaches precipitation temperature.
- Total Failure Fraction for Debris Outside the ZOI: CASA Grande uses table of total failure fractions applied to transport logic trees. Fraction of each type (fiber, paint and coatings, etc.) that passes to the pool are used to understand what is in the pool as a function of time during recirculation.
Total failure fraction multiplies total inventory of unqualified coatings.
- Chemical Bump Up Factor: Used as a multiplier on conventional head loss calculated in CASA Grande. Multiplier is applied if thin bed is formed and pool temperature is at or below precipitation temperature.
Step 3: Select Inputs of Interest
- Fiber Penetration Function: Fraction of fiber that bypasses the ECCS sump screen as a function of the amount of fiber on the screen.
- Size of ZOI: ZOI defined as direct function (multiplier) of break size and nominal pipe diameter; e.g., for NUKON fiber, ZOI is 17 times break diameter. ZOI is spherical unless break is not DEGB, in which case it is hemispherical. Truncated by any concrete walls within the ZOI.
- Time to Turn Off One Spray Pump: If three spray pumps start, then by procedure one is secured. Time to secure the pump is governed by operator acting on the conditional action step in procedure.
Step 3: Select Inputs of Interest
- Time to Hot Leg Injection: Similar to the spray pump turn off time, the time to switch one or more trains to hot leg injection operation is governed by procedure.
- Strainer Buckling Limit: Limit is the differential pressure across ECCS strainer at which strainer is assumed to fail mechanically. This limit is based on engineering calculations that incorporate safety factor.
- Water Volume in the Pool: Depending on break size, amount of water in pool, as opposed to amount in RCS and other areas in containment, varies. Smaller breaks tend to result in less pool volume than larger breaks.
Step 3: Select Inputs of Interest
- Debris Densities: Depends on amount and type of debris that arrives in pool. These densities are used in head loss correlations to calculate, for example, debris volume.
- Time Dependent Temperature Profiles: Temperature of water in sump affects air release and vaporization during recirculation. Time-dependent temperature profile comes from coupled RELAP5-3D and MELCOR simulations depending on break size.
- Spray Transport Fractions for Debris Outside ZOI: CASA Grande uses a table of failure fractions applied to transport logic trees. Fractions of each type of debris that passes to pool are used to understand what is in pool as function of time during recirculation. The spray failure fraction is fraction of failed coatings that wash to pool during spray operation.
Step 4: Nominal Values and Ranges for Inputs
&#! #! ' 0 ' 1 ' 2 ' 3 )$) '7 B8 @A4D E4AD AD D?
%'%$ ' # )7 58 F4D C @D D? ' ('$(&%') $( ( %. )'%",# $ %%" ( 6@?K L@?K # " ' & )*%$#&7%8 @C?% @E?% %)" ",'K ,)( @??K G?K # ",#&6&)%' ( LD?K ' $)'*%$,$*%$ ( 1 ( 6BBK ,'$%@&'0 ,#&7# $48 A? @CC? %) $ *%$7# $48 BCD CD? )' $',!" $ # )7 4A 8 H4BD H4E ' ($( )07"#5 B8 ( LADK #&'),' '%"(7%8 ( 6DK &'0'$(&%')K ,)( EK @AK
Step 5: Estimate Outputs Under Nominal Values of Inputs . "$'#* "&! ) # 548 548 548 8
!$ -# ( # ! #
? (" $ %$ @4G@F6?G @4H@C6?H @4EAE6?G A4??H6?G @?4DBK
Step 6: One-Way Sensitivity Analysis . ) # /,54- ' 8
"$'#* "&!
!$
? (" $ @4G@F6?G @?4DBK @ )$) ' %.7E4AD B8 '( @4H?D6?G @?4@AK A )$) ' 7AD B8 $'( @4EEH6?G @?4E@K B )$) ''0 7D? B8 $'( B4BHC6?G CA4EBK C %'%$," # ) %.7C4? 58 $'( @4EH?6?E EF4FHK D %'%$," # )'0 7D? 58 '( @4B?G6?G @?4G?K E %'%$," # ) 7@D 58 '( @4BAH6?G @?4EDK F ' ('$(&%') $( $'( F4GHE6?G AG4D?K G ' ('$(&%') $( %. '( @4AC@6?G @A4?BK H # "#& $'( @4H?D6?G @?4@FK @? %)" ",'K ,)( %.7G?K8 '( @4FF?6?G @?4E@K @@ ,#&)%' $'( A4AGF6?G G4GDK @A $)'*%$ %.$-"%& $'( @4DDA6?F @?4HBK @B 1#"" '( E4FHD6?H @A4@GK @C ,'$ @&'0 %$ ' '( @4DEH6?G @@4ABK @D %) $ *%$ %$ ' $'( @4HEA6?G H4HEK @E )' $' # ) ' '( @4EBH6?G @?4HHK @F )'%",# %. $'( A4??@6?G @?4@BK @G )'%",# '( @4EDD6?G @?4FBK @H ' ($( )0 $'( A4DEF6?G H4@FK A? #&'),' '%"( %. $'( @4HEB6?G @?4@CK A@ &'0'$(&%')K ,)( 7@AK8 $'( @4FHG6?G @?4EDK
Step 6: One-Way Sensitivity Analysis: Sensitivity Plots & Tornado Diagrams
"$!" $ - %
& "" % '!$&!$ "%&" % '!$ " ! $ &$ '%
6.76 7.66 76.66 766.66 7,666.66
"$"!' &1:.6/0;6/2
!&$& "! ")!( "#
$ %$!%#"$& !% $% $ &$ '%
+ !$% $ &$ '%
!$ "%$ #& "!%
&!& $1<.8;&390;6&392
$ %!% &*
' #&"$
&$" '
'$! 7#$* "!$
&$ !$ & $
#$&'$ $" % ")
"& ! & "! "!$
#
"& '$@ '&% ")1=6@2
#$*$!%#"$&@ '&% 178@2
Step 6: One-Way Sensitivity Analysis: Sensitivity Plots & Tornado Diagrams
$### #
,$ # (
$## #
$# #
$ #
# $
% # % ' & # & ' ' # ' ' ( # ( ' ) # ) ' * # * ' + #
}}