Topics in Model Validation and Uncertainty Quantification, Volume 5

River Rapids Conference Proceedings of the Society for Experimental Mechanics Series Topics in Model Validation and Uncertainty Quantification, Volume 5 Todd Simmermacher Scott Cogan Babak Moaveni Costas Papadimitriou Proceedings of the 31st IMAC, A Conference on Structural Dynamics, 2013 River Publishers

Conference Proceedings of the Society for Experimental Mechanics Series Series Editor TomProulx Society for Experimental Mechanics, Inc., Bethel, CT, USA

River Publishers Todd Simmermacher • Scott Cogan • Babak Moaveni Costas Papadimitriou Editors Topics in Model Validation and Uncertainty Quantification, Volume 5 Proceedings of the 31st IMAC, A Conference on Structural Dynamics, 2013

Published, sold and distributed by: River Publishers Broagervej 10 9260 Gistrup Denmark www.riverpublishers.com ISBN 978-87-7004-877-4 (eBook) Conference Proceedings of the Society for Experimental Mechanics An imprint of River Publishers © The Society for Experimental Mechanics, Inc. 2013 This work is subject to copyright. All rights are solely and exclusively licensed by the Publisher, whether the whole or part of the material is concerned, specifically the rights of translation, reprinting, reuse of illustrations, recitation, broadcasting, or reproduction in any other physical way, and transmission or information storage and retrieval, electronic adaptation, computer software, or by similar or dissimilar methodology now known or hereafter developed. The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use. The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Preface Topics in Model Validation and Uncertainty Quantification represents one of seven volumes of technical papers presented at the 31st IMAC, A Conference and Exposition on Structural Dynamics, 2013 organized by the Society for Experimental Mechanics, and held in Garden Grove, California February 11–14, 2013. The full proceedings also include volumes on Nonlinear Dynamics; Experimental Dynamics Substructuring; Dynamics of Bridges; Dynamics of Civil Structures; Special Topics in Structural Dynamics; and, Modal Analysis. Each collection presents early findings from experimental and computational investigations on an important area within Structural Dynamics. Model Validation and Uncertainty Quantification is one of these areas. Model Verification and Validation (V&V), refers to technology developed to assess the accuracy of numerical predictions and quantify the uncertainty and confidence with which predictions are made. V&V includes activities such as code and calculation verification, the design of validation experiments, test-analysis correlation, uncertainty propagation and quantification for both models and experiments, and finite element model updating. The organizers would like to thank the authors, presenters, session organizers, and session chairs for their participation in this track. Albuquerque, NM, USA Todd Simmermacher Besanc¸on, France Scott Cogan Medford, MA, USA Babak Moaveni Volos, Greece Costas Papadimitriou v

Contents 1 Optimal Inequalities to Bound a Performance Probability......................................................... 1 Franc¸ois M. Hemez and Christopher J. Stull 2 Remaining Fatigue Life Predictions Considering Load and Model Parameters Uncertainty .................. 17 Maurizio Gobbato, Joel P. Conte, and John B. Kosmatka 3 Fast Computing Techniques for Bayesian Uncertainty Quantification in Structural Dynamics ............... 25 Costas Papadimitriou and Dimitra-Christina Papadioti 4 Bayesian Uncertainty Quantification and Propagation in Nonlinear Structural Dynamics .................... 33 Dimitrios Giagopoulos, Dimitra-Christina Papadioti, Costas Papadimitriou, and Sotirios Natsiavas 5 Probabilistic Damage Identification of the Dowling Hall Footbridge Using Bayesian FE Model Updating... 43 Iman Behmanesh and Babak Moaveni 6 Considering Wave Passage Effects in Blind Identification of Long-Span Bridges ............................... 53 S. Farid Ghahari, M. Ali Ghannad, James Norman, Adam Crewe, Fariba Abazarsa, and Ertugrul Taciroglu 7 Quantification of Parametric Model Uncertainties in Finite Element Model Updating Problem via Fuzzy Numbers ...................................................................................................... 67 Yildirim Serhat Erdogan, Mustafa Gul, F. Necati Catbas, and Pelin Gundes Bakir 8 Quantifying Maximum Achievable Accuracy of Identified Modal Parameters from Noise Contaminated Free Vibration Data.................................................................................... 75 Eric M. Hernandez 9 Using P-Box and PiFE to Express Uncertainty in Model Updating................................................ 81 Ramin Madarshahian, Juan M. Caicedo, and Boris A. Za´rate 10 Robust Model Calibration with Load Uncertainties................................................................. 89 D. Pereiro, S. Cogan, E. Sadoulet-Reboul, and F. Martinez 11 Simulating the Dynamics of the CX-100 Wind Turbine Blade: Model Selection Using a Robustness Criterion.................................................................................................. 99 Kendra L. Van Buren, Sez Atamturktur, and Franc¸ois M. Hemez 12 Defining Coverage of a Domain Using a Modified Nearest-Neighbor Metric..................................... 113 Matthew C. Egeberg, Sez Atamturktur, and Franc¸ois M. Hemez 13 Orthogonality for Modal Vector Correlation: The Effects of Removing Degrees-of-Freedom................. 123 Michael L. Mains 14 CAE Model Correlation Metrics for Automotive Noise and Vibration Analysis ................................. 135 Qijun Zhang, Shawn Hui, and Kurt Schneider 15 Damage Localization Using a Statistical Test on Residuals from the SDDLV Approach........................ 143 L. Marin, M. Do¨hler, D. Bernal, and L. Mevel vii

viii Contents 16 Robust Tolerance Design in Structural Dynamics ................................................................... 153 Chaoping Zang, Jun Yang, and M.I. Friswell 17 Uncertainty Propagation in Floating Raft System by FRF-Based Substructuring Method for Elastic Coupling ..................................................................................................... 165 Huang Xiuchang, Hua Hongxing, Chen Feng, and Xu Shiyin 18 Crossing and Veering Phenomena in Crank Mechanism Dynamics ............................................... 175 Elvio Bonisoli, Gabriele Marcuccio, and Carlo Rosso 19 Validating Low-Level Footfall-Induced Vibration Predictions in Steel and Concrete Structures.............. 189 Michael J. Wesolowsky, Julia M. Graham, J. Shayne Love, Jon K. Galsworthy, and John C. Swallow 20 Finite Element Model Updating of an Assembled Aero-Engine Casing........................................... 199 Chaoping Zang, Shuangchao Ma, and M.I. Friswell 21 Experimental Modal Analysis and Modelling of an Agricultural Tire ............................................ 213 F. Braghin, F. Cheli, S. Melzi, S. Negrini, and E. Sabbioni 22 International Space Station Modal Correlation Analysis ........................................................... 221 Kristin Fitzpatrick, Michael Grygier, Michael Laible, and Sujatha Sugavanam 23 Numerical Modeling of Vibration Induced Atomization of Liquids ............................................... 243 Jesi Ehrhorn and William Semke 24 Dynamical Modeling and Verification of Underwater Acoustic System........................................... 255 Ahmet Levent Avs¸ar, ˙Istek Tatar, and Cihangir Duran

Chapter1 Optimal Inequalities to Bound a Performance Probability Franc¸ois M. Hemez and Christopher J. Stull Abstract A challenging problem encountered in engineering applications is the estimation of a probability-of-failure based on incomplete knowledge of the sources of uncertainty and/or limited sampling. Theories formulated to derive upper probability bounds offer an attractive alternative because first, they avoid postulating the probability laws that are often unknown and second, they substitute numerical optimization for statistical sampling. A critical assessment of one such technique is presented. It derives upper probability bounds from the McDiarmid concentration-of-measure theory, which postulates that fluctuations of a function are more-or-less concentrated about its mean value. Two applications of this theory are presented. The first application analyzes a “toy” polynomial function defined in two dimensions. The upper bounds of probability are calculated and compared to sampling-based estimates of the true-but-unknown probabilities. For this function, the upper bounds obtained are too broad to be useful. These results are confirmed by conducting a similar analysis on a real engineering system, where upper bounds of probability associated with resonant frequencies of a structural system are estimated. A high-fidelity finite element model, previously validated using vibration measurements, is used to predict the frequencies. In this application, the uncertainty is introduced by way of material properties and the effective preload of a beam-to-column connection, modeled explicitly. These applications suggest that the theory not only leads to upper bounds that are inefficient but that can also be sub-optimal if their numerical estimation is based on too few model runs. It is concluded that this particular theory, while mathematically attractive, may not be well suited for engineering applications. Keywords Concentration-of-measure • Probability-of-failure • McDiarmid diameter • Parametric uncertainty • Model validation (Approved for unlimited, public release on September 15, 2012, LA-UR-12-24769, Unclassified.) 1.1 Introduction An analysis of system reliability estimates a probability-of-failure, denoted as PF, where “failure” is defined as a condition where the “demand” exceeds the “capacity” that an engineered system is capable of providing. Demand, for example, refers to a peak stress resulting from an applied load, while capacity refers to a material yield stress criterion. Another application would be to require that a performance metric, that defines the capacity C, meets a user requirement, or demand D. The probability-of-failure can be written in a generic sense as: PF =Prob[D≥C], (1.1) F.M. Hemez ( ) Technical Staff Member, Los Alamos National Laboratory, X-Theoretical Design Division (XTD-3), PO Box 1663, Mail Stop T087, Los Alamos, NM 87545, USA e-mail: hemez@lanl.gov C.J. Stull Technical Staff Member, Los Alamos National Laboratory, Applied Engineering and Technology Division (AET-6), PO Box 1663, Mail Stop P915, Los Alamos, NM 87545, USA e-mail: stull@lanl.gov T. Simmermacher et al. (eds.), Topics in Model Validation and Uncertainty Quantification, Volume 5, Conference Proceedings of the Society for Experimental Mechanics Series 41, DOI 10.1007/978-1-4614-6564-5 1, © The Society for Experimental Mechanics, Inc. 2013 1

2 F.M. Hemez and C.J. Stull Probability Density Function Performance Demand, D Capacity, C Region that contributes to the failure probability, PF Fig. 1.1 Generic definition of a probability-of-failure PF where the symbol Prob[•] denotes a probability function. The probability-of-failure PF is shown graphically in Fig. 1.1 as the shaded area, where the Probability Density Function (PDF) of encountering a demand that exceeds capacity is non-zero. The failure probability integrates the area under the curve and the reliability R of the engineered system is, by definition, equal to: R=1−PF. (1.2) Equations (1.1) and (1.2) indicate that the probability laws of the random variables (C; D) must be known analytically or, at a minimum, can be sampled a sufficient number of times, in order to estimate PF with reasonable accuracy. In many applications, however, this is not possible. These include cases where the datasets available to estimate the pair (C; D), and their uncertainties, are too sparse to be useful for numerical integration. Another commonly encountered situation is one where the capacity C is estimated numerically by analyzing the output of a numerical simulation. The large computational expense of running the model may prevent the accumulation of enough samples to yield satisfactorily converged statistics. Also, the high dimensionality of an input space may make it difficult to numerically integrate the probability-of-failure (1.1). This report evaluates another approach to bypass the challenge of sampling a high-dimensional space. This alternative approach searches for an upper bound of the failure probability. Put forth by the Predictive Science Academic Alliance Program (PSASP) Center at Caltech in Pasadena, California, it is a technique that relies on concentration-of-measure inequalities [1]. The method derives an upper bound UP of the failure probability: PF ≤UP, (1.3) without having to assume probability laws for the demand and capacity terms of Eq. (1.1). Such an upper bound can be derived from a variety of theories. The technique of Ref. [1] is based on the concentration-of-measure theory that postulates that “fluctuations” of a function are more-or-less concentrated about the mean value of the function. This assumption is taken advantage of to derive the upper bound of probability UP. The technique is evaluated because of the great potential that it offers. As alluded to, the upper bound is by construction, valid irrespective of the probability laws of the variables (D; C). This is an attractive property in cases where the datasets analyzed are so sparse that probabilities cannot be inferred reliably. It also is useful in situations where sampling the function is computationally expensive. The examples discussed herein evaluate the upper bounds when applied to nonlinear functions that are defined with two or three variables and fluctuate significantly within their design spaces. Nonlinearity, moderate dimensionality, and large fluctuations characterize the class of function properties of interest for engineering applications. The material modeling discussed in Ref. [2], for example, involves 7 parameters while the simulations of Refs. [3–5] are defined with 20 variables, approximately. Higher-dimensional problems are, at this time, considered to be out-of-reach for practical applications. Unfortunately, the upper probability bounds derived using concentration-of-measure inequalities do not live up to their promises, at least not for the examples considered in this work. The main two observations are that: first, the upper bounds are too broad to be useful for sensible decision-making; and second, the computational cost of their estimation is prohibitive even for the relatively low-dimensional problems considered. Furthermore, it is observed that the true-but-unknown failure probability PF can exceed the upper bound UP if its estimation relies on too few samples or iterations, which can lead to a dangerous, false sense of confidence. After briefly overviewing the derivation of upper probability bounds in Sect. 1.2, two examples are presented. Section 1.3 discusses a “toy” problem where the upper bounds are estimated for the Rosenbrock function in two dimensions: a polynomial whose numerical optimization is known to be challenging. The discussion includes a statistical convergence study of UP that underscores the danger of “under-sampling.” The second application, discussed in Sect. 1.4, analyzes a

1 Optimal Inequalities to Bound a Performance Probability 3 more realistic engineering problem that involves a high-fidelity finite element simulation of structural vibration response. Section 1.5 concludes the discussion and offers recommendations. 1.2 Concentration-of-Measure Inequalities of Probability Bounds This section gives a brief summary of the derivations presented in Refs. [1,6] to establish the upper bound of a probability-offailure. The approach is based on concentration-of-measure inequalities [7]. We do not attempt to re-derive these results and focus instead on the equations needed to understand the numerical implementations and applications of Sects. 1.3 and 1.4. 1.2.1 Upper Probability Bounds Based on Concentration-of-Measure Inequalities The proposal of Ref. [1] for uncertainty quantification and system certification is that the fluctuation of a function is more-orless concentrated about its mean value. The extent of this concentration is controlled by a “diameter” quantity DF that canbe estimated experimentally or numerically. This concept of concentration is exploited to bound the probability that a function exceeds its expected (or mean) behavior. Here, the function represents a complex engineered system whose performance we wish to evaluate. The terminology “deviation from the mean,” likewise, indicates a variation away from the nominal performance of the system. The probability of exceeding a given fluctuation represents the failure probability of the system. The main result that is illustrated with several applications is that the probability of observing a fluctuation of a given magnitude away from the mean value is bounded according to: Prob[F(X)−μ≥ δ] ≤e− 2 δ D F 2, (1.4) where F(•) is a real-valued function that depends on N input variables X=(X1;X2;...XN) and δis the level of fluctuation away from the mean value, denoted by μ, of the function: μ≡E[F]= X∈([−1;+1]) N F(X)· dX. (1.5) Equation (1.5) implies a multi-dimensional integration in the space ℜN. Without loss of generality it is assumed that the variables Xk are scaled in the interval −1≤Xk ≤+1. The symbol DF of Eq. (1.4) denotes the McDiarmid diameter of the function that controls the extent to which fluctuations are concentrated (or not) about the mean value. To derive an upper bound of probability at any level of fluctuation δ, it suffices to know the diameter DF that is formally defined as: D2 F = ∑ 1≤k≤N D2 k, (1.6) whereDk is a partial diameter that defines the influence on function values of the kth input variable Xk: Dk = sup {X,X∗∈U(k)} (F(X)−F(X∗)). (1.7) Equation (1.7) seeks the maximum range of function F(•) over a specific space U(k). This space, within which the optimization is performed, is the subspace of the hyper-cube domain ([−1;+1])N where the kth input variable Xk is constrained to be constant: U(k) = X∈([−1;+1]) Nsuch thatXk =constant . (1.8) Note that keeping the kth input variable constant and equal to a known value leads to a partial diameter Dk that becomes conditioned on Xk. This is not what is meant in Eqs. (1.7) and (1.8). Instead, the optimization is defined in the N-dimensional space ([−1;+1])N while enforcing the constraint that Xk remains constant and equal to an unknown value. This unknown value ofXk must also be optimized while solving for Eq. (1.7).

4 F.M. Hemez and C.J. Stull 1.2.2 Numerical Implementation of the Technique Often, an engineered system is designed such that its performance does not deviate “too much” from the nominal (or mean) target performance. “Failure” can therefore be conceptually defined as exceeding a fluctuation of response that would be detrimental to the required performance. The corresponding probability-of-failure is: PF =Prob[F(X)−μ≥ δ], (1.9) and the reliability of the engineered system is defined, as before, by Eq. (1.2). Equations (1.4) through (1.6) define a rather simple procedure to derive an upper bound of reliability. It rests on the ability to evaluate, either experimentally or numerically, the McDiarmid diameter. Equation (1.6) indicates that DF accumulates the partial diameters Dk in quadrature; the simplicity of this definition, however, is deceiving. Estimating a partial diameter Dk can be viewed, in fact, as a two-level, nested optimization. The innermost optimization optimizes the kth variable while the other variables are kept constant and equal to fixed values, that is, X1 =λ1, X2 =λ2,...Xk −1 =λk−1 , Xk+1 =λk+1,...XN =λN. The variable Xk is optimized to search for the maximum range F(Max) k (λ(−k))− F(Min) k (λ(−k)) shown below, where the symbol λ( −k) denotes the vector of ℜN−1 defined as λ(−k) =(λ1 ; λ2;...λk −1 ;λk+1;...;λN). (The subscript “-k” is used to indicate that all variables are included, except the kth one.) This first step is embedded in an outer optimization that solves for the partial diameter Dk by optimizing the (N−1) variables of vector λ(−k) . These two steps must be repeated for each variable. The numerical implementation of the nested optimizations (1.7) and (1.8) can be written as: Dk = max {λ(−k)∈([− 1;+1]) N−1} F(Max) k (λ( −k))− F(Min) k (λ( −k)), (1.10) where F(Min) k (λ( −k))= min {−1≤Xk≤+1} F(λ1; λ2;··· λk −1 ;Xk; λk+1;··· λN), (1.11) and F(Max) k (λ( −k))= max {−1≤Xk≤+1} F(λ1; λ2;··· λk −1 ;Xk; λk+1;··· λN). (1.12) It can be verified that substituting the minimum (1.11) and maximum (1.12) in Eq. (1.10) is equivalent to searching for the maximum range, written as “sup(F(X)−F(X∗))” in Eq. (1.7). Equations (1.10) through (1.12) reveal that the calculation of a partial diameter Dk requires the resolution of three optimization problems over different subspaces of ([−1;+1]) N. This is a rather expensive proposition since the resolution of these optimizations must be repeated for every variable X1, X2,...,XN of the problem. Three strategies are implemented to estimate the McDiarmid diameter. The first two algorithms, applied in Sect. 1.3 to the two-dimensional Rosenbrock polynomial function and in Sect. 1.4 to the three-dimensional structural vibration problem, rely on the formulation of Eqs. (1.10) through (1.12). In the first strategy, the functions F(Min) k (•) and F (Max) k (•) are explored for several values of λ( −k) by solving the minimization and maximization problems of Eqs. (1.11) and (1.12). After enough “samples” have been obtained, multidimensional emulators are best-fitted to the pairs (λ( −k) ; F(Max) k (λ(−k))) and (λ( −k); F(Min) k (λ(−k))). The emulators are defined as polynomial models in subspace ([−1;+1])N−1. Their orders are kept as low as possible, provided that the goodness-of-fit is acceptable, to avoid over-fitting the training data. These fast-running emulators are used to solve the optimization of Eq. (1.10). The procedure is then repeated for every variable Xk of the problem. The second strategy adopts a nested implementation of Genetic Algorithms (GAs) to optimize Eq. (1.10). The GA-based implementation provided by the Global Optimization Toolbox in MATLABTM [8] is employed, where it is noted that the default settings with respect to selection, crossover, and mutation functions are used for the analyses discussed herein. As GAs are well-accepted in the numerical optimization community, only the details relevant to the nested implementation are provided. First, the “outer” GA initiates a population of individuals. Each of these individuals is then passed to an “inner” GA that searches the space associated with the kth parameter to determine the maximum range of the function, while keeping the λ(−k) variables fixed. That is, each execution of the inner GA seeks the maximum range for a fixed subset of λ( −k) variables, instead solving Eqs. (1.11) and (1.12) separately. After all individuals from the outer GA are evaluated in this manner, they are ranked, selected, etc., and the procedure repeats for each population of individuals until termination of the outer GA.

1 Optimal Inequalities to Bound a Performance Probability 5 The third approach implemented for comparison with the first two consists of relying on a large number of Monte Carlo samples, then, searching for the minimum and maximum values of functions F(Min) k (•) andF (Max) k (•). The implementation of this sampling strategy is as follows. First, a predefined number of samples M, are generated from the space ([−1;+1])N. This produces M vectors of randomly sampled variables, denoted as the “unprimed” samples. Next, for each of the M vectors, only the kth variable is re-sampled; these new vectors are denoted as the “primed” samples. The final step is then to evaluate the unprimed and primed samples and store the absolute value of the difference between the mth unprimed and the mth primed evaluation (where m=1,2,...,M). The maximum difference found among the M samples is then designated as the partial diameter Dk. In this way, the two-step optimization is being condensed into a single step, whereby the λ( −k) variables are held constant for two random (i.e. unknown) values of the kth variable. The quality of the estimation rests entirely on whether or not the two random values of the kth variable generate the minimumandmaximum values of functions F(Min) k (•) and F(Max) k (•), for the “correct” values of the λ( −k) variables. Lastly, in order to assess the convergence properties of the sampling strategy, the above procedure is repeated ten times. 1.3 Application to the 2D Rosenbrock Function The first application of the concentration-of-measure inequality-based reliability is to a two-dimensional polynomial function. It represents a “toy” problem that does not reflect the practical difficulties of real engineered applications. This problem, nevertheless, offers the advantage of simplicity. It also allows for the assessment of convergence properties, discussed in this section. 1.3.1 The Two-Dimensional Rosenbrock Function The Rosenbrock function is a simple, analytic polynomial whose numerical evaluation is trivial. The results discussed are obtained with an implementation developed with the software MATLABTM [8]. After estimating the McDiarmid diameter described by Eq. (1.6), the upper bound in Eq. (1.4) is calculated where “failure” is defined as the probability that the Rosenbrock function deviates from the mean value defined by Eq. (1.5) by more than a given fluctuation δ. The results are presented as curves that show the upper probability bounds as a function of δ. The Rosenbrock function in the N-dimensional space (X1; X2;...XN) is defined analytically as: Y = ∑ k=1···(N−1) (1−Xk) 2 +100· Xk+1−X 2 k 2 , (1.13) where each variable Xk is scaled in [−1; +1]. In two dimensions, the function becomes: Y =(1−X1) 2 +100· X2−X 2 1 2 . (1.14) The Rosenbrock function is a benchmark test commonly encountered in numerical optimization because it varies rapidly and defines “valleys” [9]. This makes it challenging to search for a global minimum using gradient-based optimization algorithms. Another reason for selecting this function is the existence of previous work documented in Refs. [10, 11]. Reference [12] presents results obtained with a variant defined in 17 dimensions. Figure 1.2 illustrates the Rosenbrock function in two dimensions. The formation of a “valley,” that makes numerical optimization difficult, can be observed to be centered at the coordinates (X1; X2)=(0;0). The valley bifurcates into two branches: one of them progresses towards the corner (X1;X2)=(+1;−1)while the other branch extends towards(X1;X2)= (+1;+1). It can be verified that the global minimum of the Rosenbrock function is reached when Xk =+1 for all variables, whether in 2D or for an arbitrary number of dimensions, which is the only point where the function is equal to zero. In two dimensions, the maximum of the function over ([−1;+1])2 is equal to 404. Another theoretical result [7] combines the minimum and maximum values to bound the McDiarmid diameter as: 1 √N≤ DF max {X∈([−1;+1]) N} F(X)− min {X∈([−1;+1]) N} F(X) ≤ √N. (1.15)

6 F.M. Hemez and C.J. Stull Fig. 1.2 Illustration of the Rosenbrock function in 2D space The numerical application of Eq. (1.15) indicates that the McDiarmid diameter is bounded in the interval 285.67≤DF ≤ 571.34. These bounds, although they are somewhat “loose,” can be used to assess the quality of estimations obtained with different numerical implementations. The next statistics useful for the analysis are the mean and standard deviation. In 2D, the mean value ( μY) and standard deviation ( σY) are obtained analytically by integrating the function over the hyper-cube space ([−1;+1])2. The intermediate derivations for these statistics are omitted for brevity and instead given as: μY = 1 Volume [−1;+1] 2 · +1 −1 +1 −1 (1−X1) 2 +100· X2−X 2 1 2 dX1dX2 ≈54.667, (1.16) and σY = 1 Volume [−1;+1] 2 · +1 −1 +1 −1 (1−X1) 2 +100· X2−X2 1 2 −μY 2 dX1dX2 ≈65.433. (1.17) It is noted that the standard deviation of the 2D Rosenbrock function is large relative to its mean value, due to the presence of corners where values of the polynomial become large (see Fig. 1.2). This is a property encountered in other engineering applications but is not important to justify the method, as indicated in Ref. [12] where it does not hold in the case of the 17D Rosenbrock variant. The mean and standard deviation values, μY and σY, are used to guide the selection of appropriate function fluctuations δin the analysis that follows. 1.3.2 Estimation of the McDiarmid Diameter and Upper Probability Bounds The analysis begins by estimating the true-but-unknown failure probabilities using statistical sampling. Collecting this information is, here, possible because the Rosenbrock function can be evaluated numerically in a fraction of second. These highly accurate estimations of probability are used for comparison with the McDiarmid-based upper bounds of Eq. (1.4). Figure 1.3 illustrates the empirical histogram of Rosenbrock function values. It is obtained from ten million (10 M) uniform and uncorrelated samples of variables (X1; X2) in hyper-cube ([−1;+1]) 2. Analyzing 10 M samples suffices to yield statistically converged estimates, an assertion that is verified by repeating the estimation ten times and observing that the statistical estimates do not change in their fourth significant digit from one trial to the next. The resulting statistics are

1 Optimal Inequalities to Bound a Performance Probability 7 Fig. 1.3 Histogram of values sampled for the 2D Rosenbrock function Table 1.1 Statistics of the 2D Rosenbrock function estimated with 10 M samples Definition Symbol Estimate “Exact” value Minimum value Y(Min) 0.00 0.00 Maximum value Y(Max) 404.00 404.00 Total range Y(Max) −Y(Min) 404.00 404.00 Mean value μY 54.67 54.67 Standard deviation value σY 65.43 65.43 listed in Table 1.1, where they are compared to the exact values presented in Sect. 1.3.1. The “exact” failure probabilities are estimated from these samples in a frequentist manner, simply by counting how many values exceed a given level of fluctuation δrelative to the mean value μY. Next, the conditional functions F(Min) k (•) andF (Max) k (•) of Eqs. (1.11) and (1.12) are estimated numerically. Because the function is defined in two dimensions, the symbol λ( −k) that denotes the combination of variables kept constant during the optimization reduces to a single variable. The conditional functions F(Min) 1 (•) and F (Max) 1 (•) represent the minimum and maximum of the Rosenbrock function, respectively, as variable X2 is kept constant and fixed to a known value. Likewise, F(Min) 2 (•) and F (Max) 2 (•) are the minimum and maximum as variable X1 is kept constant. A total of 50,000 equally-spaced samples are selected in the interval [−1; +1] to generate these constant values “X2 =constant” (for F (Min) 1 and F (Max) 1 ) or “X1 =constant” (for F (Min) 2 andF (Max) 2 ). Each time that one of these fixed values is selected for one of the variables, the other variable of the function is optimized to search for F(Min) k andF (Max) k . It is emphasized that these are simple, one-dimensional constrained optimization problems. They must, however, be repeated twice to search for the minimum and maximum values for each “sample” considered (which here, is 50,000 evaluations). Once F(Min) k and F (Max) k have been evaluated numerically a sufficient number of times over the design space, the partial function range Gk(•) is defined as: Gk(λ( −k))=F(Max) k (λ( −k))− F(Min) k (λ( −k)). (1.18) Because the partial function range Gk(•) depends on a single variable, its functional dependency on X( −k) can be suggested graphically. The run numbers are sorted in ascending magnitudes of X( −k) and the resulting values are shown in Fig. 1.4a, b forG1(•) andG2(•), respectively. Figure 1.4a shows that the partial function range G1(X2) is smoothly varying, with the exception of a sharp discontinuity at X2 =0.19. Figure 1.4b suggests a quadratic functional form for G2(X1). Going back to Eq. (1.10), it can be inferred that the partial diameter Dk is simply the maximum of function Gk(•), which can be read directly from the figure. Of course, higher-dimensional problems necessitate a numerical optimization to estimate Dk from the function Gk(•).

8 F.M. Hemez and C.J. Stull Fig. 1.4 Sorted training data and polynomial emulators. (a) Left: partial function range G1(X2) versusX2. (b) Right: partial function range G2(X1) versusX1 Table 1.2 “Exact” and approximate values of the McDiarmid partial diameters, Dk McDiarmid partial diameter, Dk (no unit) Input variable, Xk “Exact” value of the diameter Optimization-based estimation D1 303.00 123.79 D2 400.00 400.00 DF 501.80 418.72 The next step is to develop a low-order polynomial emulator from the training data accumulated for each partial function range Gk(•). This step is, here, not needed since the maximum values Dk can be read directly from Fig. 1.4a, b; the polynomial emulators are nevertheless developed for consistency when solving higher dimensional problems as discussed in [12]. Fourth-order polynomials are best-fitted to the 50,000 evaluations of functions G1(•) andG2(•). The goodness-of-fit obtained for G1(•) is 8.60% Mean Square Error (MSE). Given that the function exhibits a discontinuity, this level of accuracy is deemed appropriate. The goodness-of-fit of the fourth-order polynomial developed to emulate G2(•) is 5.29×10− 7% MSE. Figure 1.4a, b visually compare the polynomial emulators to the training data. Clearly, the discontinuity of function G1(•) in Fig. 1.4a cannot be approximated by a polynomial, irrespective of its order. What matters most, however, is that the maximum of the data be represented appropriately because it is what the optimization is after. The surrogate models shown in Fig. 1.4a, b are optimized to search for the maxima, which are the partial diameters Dk. Results of the optimization are tabulated in Table 1.2 where they are compared to the exact solutions. It is observed that the first partial diameter D1 is severely under-estimated; the second optimization converges to the correct value for D2. The McDiarmid diameter DF is then calculated from Eq. (1.6) as: DF = D2 1+D2 2 = (123.79) 2 +(400.00) 2 =418.72. (1.19) This estimate complies with the interval 285.67≤DF ≤571.34 of Eq. (1.15). Note that for the remainder of this discussion, the valueDF =418.72 is used, as opposed to the correct answer listed in Table 1.2. The upper bounds of failure probability defined in Eq. (1.4) are obtained for different levels of fluctuation δusing the mean statistic of Eq. (1.16) and McDiarmid diameter of Eq. (1.19). Figure 1.5 compares the upper bounds to the true-but-unknown probabilities estimated from 10 M uniform and uncorrelated Monte Carlo samples. Levels of fluctuation ranging from 0 to 202 are used to display the results in Fig. 1.5. A fluctuation of δ=202 represents three standard deviations away from the mean value of the Rosenbrock function. It can be observed that, irrespective of the magnitude of fluctuation δ used, the upper bounds are too conservative to be useful for sensible decision-making. This is expected near the mean value, μY =54.67, where the concentration-of-measure theory may not apply. But it should not be the case in the “tail” of the distribution. For example, the upper bound is UP =61.6% at the last level of fluctuation, δ=202, while the true probability is only PF =2.2%.

1 Optimal Inequalities to Bound a Performance Probability 9 Fig. 1.5 Upper bounds UP and estimated probabilities PF of the 2D Rosenbrock function 0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 60 70 80 90 100 Level of Fluctuation, rY Probability-of-Failure, PF=Prob[Y(X)-μY>rY] (%) Convergence of McDiarmid Upper Bound Actual 200 2000 20000 200000 2e+06 Optimal Fig. 1.6 Sampling-based upper probability bounds UP of the 2D Rosenbrock function 1.3.3 Convergence Behavior of Sampling-Based Estimates As mentioned at the end of Sect. 1.2.2, two other strategies are available to estimate the partial diameters Dk. The second approach is based on the nested GA optimization described previously and the third approach is based on sampling the hyper-space ([−1;+1])N, assuming uniform and uncorrelated probability distributions for the variables Xk, to search for the minimum and maximum functions F(Min) k (•) and F (Max) k (•). A comparison with the optimization-based results reported Fig. 1.5 is briefly discussed next. Figure 1.6 compares the upper bounds UP of the GA optimization-based solution (dotted line) to five sampling-based estimations obtained with different numbers of Monte Carlo runs (dashed lines). Including the replicates, the sizes of these sampling-based simulations (i.e. 10×M) vary from M=200 toM=2 million samples, where it is noted that Fig. 1.6 reflects

10 F.M. Hemez and C.J. Stull 0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 60 70 80 90 100 Comparison of Probability-of-Failure and McDiarmid Upper Bound Level of Fluctuation, rY Probability-of-Failure, PF=Prob[Y(X)-μY>rY] (%) Prob[Y(X)-μY>rY] exp(−2*(rY/DF,min) 2) exp(−2*(rY/DF,avg) 2) exp(−2*(rY/DF,max) 2) Fig. 1.7 Statistics of sampling-based upper probability boundsUP obtained with 200 runs only the “best” sampling-based estimation of the upper bound UP (i.e. the maximum McDiarmid diameter) realized from the ten replicates. The “true” failure probabilities PF, estimated via Monte Carlo sampling, are also included for reference (solid, black line). Comparing Figs. 1.5 and 1.6, it is seen that both optimization-based implementations yield very nearly the same upper bounds. However, a somewhat counterintuitive result is also seen in that increasing the number of samples produces upper bounds that converge to a solution that does not correspond to either of the optimization-based solutions. This is almost entirely due to the underestimation of the first partial diameter D1, which is a common feature of both optimization-based implementations. A second observation that is somewhat more disconcerting is that the upper bounds obtained through sampling can be smaller than the values provided by the optimization-based solver. This behavior opens the door to the possibility of estimating a value UP thatwouldnot be an upper limit, therefore violating Eq. (1.3). Figure 1.7 illustrates that this dangerous situation can manifest itself when only a 200 samples are used, where it is seen that the “worst” sampling-based estimation of the upper bound (i.e. the minimum McDiarmid diameter) falls belowthe estimate of true probability-of-failure. This is an unfortunate observation as it suggests that analyzing a small number of samples may not be a viable alternative to an optimization-based solution. 1.4 Application to the High-Fidelity Model of a Three-Story Frame In this section, the upper probability bounds are applied to an engineering application. The example is the linear vibration of a three-story frame structure simulated numerically with a high-fidelity Finite Element (FE) model. The number of variables is not significantly different from the dimensionality of the Rosenbrock function discussed previously. The main difficulty of this application is to estimate the McDiarmid-based inequality, given that the FE model is computationally expensive. 1.4.1 Description of the Three-Story Frame Structure The three-story structure of Fig. 1.8a is the system employed in this study. It consists of aluminum columns and plates assembled using bolted joints. Each floor plate measures 30.5×30.5×2.5 cm, and is connected to the floors above and/or below by four aluminum columns measuring 17.7×2.5×0.6cm. Each connection comprises four bolts that act to sandwich the end of the column between the floor plate and a 2.5×2.5×1.3 cm end cap. The structure slides on rails that permit movement in one direction only, as shown in Fig. 1.8b, which is referred to as the “weak” bending direction of the columns.

1 Optimal Inequalities to Bound a Performance Probability 11 b a Fig. 1.8 High-fidelity finite element model of the three-story frame structure. (a) Left: frame structure. (b) Right: finite element discretization Table 1.3 Definition of uncertainty variables of the frame structure Variable, Xk Definition of variable X1 Elastic modulus of all Aluminum-based components X2 Mass density of all Aluminum-based components X3 Bolt “radius of influence” for end cap connections Additionally, a small column is suspended from the third floor, and a bumper mechanism is attached to the second floor. Depending upon the amplitude of the force imparted by the shaker and the gap between the suspended column and bumper, a contact nonlinearity may be introduced into the system. This capability is not exercised here. Figure 1.8b illustrates the highly refined FE model developed using the commercially available software Abaqus 6.10-1 [13]. All floors and end caps are modeled with linear, hexagonal continuum elements (Abaqus C3D8R). The bottom supports, suspended column, and bumper assembly are modeled with modified quadratic, tetrahedral continuum elements (Abaqus C3D10M). The latter element is formulated to properly account for pressures associated with contact, an important point to consider as all contact is modeled explicitly in the model. The vertical columns spanning between the base floor and the other levels are modeled with shell elements (Abaqus S4R), due to their enhanced capability to capture bending behavior, which is the dominant action experienced by those members. Figure 1.8b indicates the high geometrical fidelity with which the three-story frame structure is discretized. The mechanics are also modeled with high fidelity since all sources of contact and friction are represented explicitly. A mesh convergence study, omitted herein for brevity, provides a mesh seed of 3 mm, which is slightly less than one-eighth of an inch. The corresponding discretization, shown in Fig. 1.8b, yields predictions of the resonant frequencies that are “converged-enough” relative to the overall level of experimental variability. Details on the modeling and analysis of the frame structure are available from Ref. [14]. 1.4.2 Upper Probability Bounds of the Three-Story Frame Structure Table 1.3 presents the uncertain variables of the 3D FE model. The first two variables are straightforward in their meaning, whereas the third variable warrants explanation. It denotes the “radius of influence” associated with the bolts that connect the columns to the floors by way of the end caps. These radii are identified as red dots in Fig. 1.8b. The effect of this variable is somewhat analogous to bolt torque. Varying the value changes the rigidity of the column-to-floor connection. That is, at small values of the radius, the column-to-floor connections behave more like pinned connections (free to rotate) whereas for larger values, the beam-to-floor connections behave more like fixed connections. At very large values, the columns begin to act as rigid structural members.

12 F.M. Hemez and C.J. Stull Fig. 1.9 Goodness-of-fit of bending frequencies predicted by the surrogate models The question illustrated in this application is: given the modeling uncertainty embodied by the variables of Table 1.3, what is the probability PF that a bending frequency of the structure deviates from the mean frequency by more than, say, three Hertz? This question can be formulated as a reliability problem, where the probability of predicting a frequency is integrated for all frequency values that exceed the mean statistic by, at least, 3 Hz. This true-but-unknown value of probability PF is estimated, first, with Monte Carlo sampling. The question is then addressed using the McDiarmid-based inequality of Eq. (1.4). Comparing the upper bound UP to the true-but-unknown probability PF assesses the usefulness of this approach for decision-making. Compared to the Rosenbrock function, the 3D FE model requires a significant computational demand to calculate the lowfrequency bending modes. To address this difficulty, a design-of-experiments is populated to generate the training datasets needed to develop fast-running, polynomial emulators. These polynomial surrogates are defined as quadratic models: Y= β0+β1X1+β2X2+β3X3+β12X1X2+β13X1X3+β23X2X3+β11X2 1+β22X2 2+β33X2 3, (1.20) and yield the predictions of resonant frequency presented in Fig. 1.9. The figure plots the bending frequencies predicted by the 3D FE model as a function of those predicted by the polynomial surrogates of Eq. (1.20). The fact that the points shown are aligned with the 45◦-angle diagonal indicates that the polynomial equations emulate well the computationally expensive FE model. The high goodness-of-fit lends credence to employing these surrogates. The optimization-based and samplingbased evaluations of the McDiarmid diameter of Eqs. (1.6) and (1.10) through (1.12) are performed using the polynomial surrogates of Eq. (1.20). Figure 1.10 shows the empirical distribution of first bending frequencies predicted using Monte Carlo sampling. This histogram is obtained by analyzing 10 M samples, replicated ten times. The mean and standard deviation statistics are μY29.38Hzand σY=2.51 Hz, which gives 8.54 % variability relative to the mean. The statistics of minimum and maximum function values are equal to 22.43 and 36.46 Hz, respectively. Using Eq. (1.15), the McDiarmid diameter should fall within the bounds as: 1 √3 ≤ DF 36.46−22.43 ≤ √3, (1.21) or8.10≤DF ≤24.30. Table 1.4 indicates the partial diameters estimated from 10 M samples, replicated ten times. Combining definition (1.6) and the results listed in the table yields a McDiarmid diameter of DF =9.67 Hz, which “clears” the inequalities of Eq. (1.21). Figure 1.11 depicts the upper probability bounds UP as a function of frequency fluctuation, δ. Solutions are obtained with five sampling sizes, ranging from M=300 to M=3 million samples. The upper bound reached by relying on numerical optimization is compared to the sampling-based bounds. Also included for comparison are the true-but-unknown probabilities obtained by integrating the empirical histogram of Fig. 1.10 for each value of δ. The first observation from Fig. 1.11 is that the sampling-based estimation of upper probability bounds converges toward a solution that does not correspond with the optimization-based solution. A similar behavior was noticed for the 2D Rosenbrock function discussed in Sect. 1.3. It is confirmed that, as in the latter case, the optimization procedure fails

1 Optimal Inequalities to Bound a Performance Probability 13 22 24 26 28 30 32 34 36 38 0 1 2 3 4 5 6 Function value, F(Y) (no unit) Probability-of-Occurrence (%) Empirical Histogram of XXX Function Values Fig. 1.10 Histogram of first bending frequencies sampled for the 3D FE model 0 1 2 3 4 5 6 7 0 10 20 30 40 50 60 70 80 90 100 Convergence of McDiarmid Upper Bound Actual 300 3000 30000 300000 3e+06 Optimal Level of Fluctuation, rY Probability-of-Failure, PF=Prob[Y(X)-μY>rY] (%) Fig. 1.11 Sampling-based upper probability bounds UP of the 3D FEmodel Table 1.4 Partial diameters Dk of the 3D FE model estimated from 10 M samples Variable, Xk Variable definition McDiarmid partial diameter, Dk X1 Elastic modulus of Aluminum components 6.56 X2 Mass density of Aluminum components 6.59 X3 Bolt “radius of influence” of end connections 2.67 to produce the optimal (i.e. largest) partial diameters (as approximated in Table 1.4), thus yielding a reduced McDiarmid diameter. A second finding from Fig. 1.11, which has been consistently observed throughout this work, is the “looseness” of upper bounds UP relative to the true-but-unknown probabilities PF. These two observations, and those from the previous sections, make it difficult to justify using the McDiarmid inequalities to bound a failure probability.

14 F.M. Hemez and C.J. Stull 1.5 Conclusion A bounding theory is investigated for its ability to derive an upper limit for a probability-of-failure. The approach investigated is based on the McDiarmid concentration-of-measure theory, which postulates that fluctuations of a function are moreor-less concentrated about its mean value. This is an attractive proposition because, first, the method avoids postulating probability laws that are often unknown and, second, it substitutes efficient numerical optimization to slow-converging statistical sampling. Two examples of the application of this theory are presented. The first one involves the two-dimensional Rosenbrock function and the second example is a realistic engineering problem wherein a finite element model is employed to simulate the vibration of a three-story frame structure. Comparisons performed between the upper bounds and highly-accurate estimates of the true-but-unknown probabilities indicate that these bounds are too “loose,” or conservative, for useful decision-making. The finding applies whether fluctuations of the function are assessed near the mean value or in the tails of the statistical distributions. In addition, it is observed that estimating the McDiarmid diameter is computationally expensive, which reduces the potential benefits if the problem studied is defined in a high-dimensional space. These observations suggest that the theory not only leads to upper bounds that are inefficient, but that can also turn out to be sub-optimal if their numerical estimation is based on too few runs of the computational model. One explanation is that the McDiarmid-based inequalities yield upper bounds that are not necessarily realized by any model of the family considered. This means that there may not be any model, included in the uncertainty space defined for analysis, whose predictions reach the upper bound UP of the inequality “PF ≤UP.” The Center of the Predictive Science Academic Alliance Program at Caltech has recently advanced another theory to provide “tight” bounds of a probability-of-failure [15–17]. The main advantage of this novel theory over the approach investigated herein is that the upper bounds, while optimal, are guaranteed to be reached by a model within the uncertainty space defined. The recent work also comes with a numerical implementation that could yield significant computational savings. Future work will investigate this novel theory, and the accompanying implementation, to provide upper bounds of a probability-of-failure. Acknowledgements This work is performed under the auspices of the Verification and Validation (V&V) program for Advanced Scientific Computing (ASC) at Los Alamos National Laboratory (LANL). The first author is grateful to Frederick J. Wysocki, V&V program manager at LANL, for his continuing support. Funding for the second author was sponsored by the U.S. Department of Energy, Nuclear Energy Division, Advanced Modeling and Simulation Office (NE-71), Nuclear Energy Advanced Modeling and Simulation (NEAMS) Program, Verification, Validation and Uncertainty Quantification (VU) Program Element. The second author is sincerely grateful to Brian J. Williams (CCS-6) for his continued support. LANL is operated by the Los Alamos National Security, LLC for the National Nuclear Security Administration of the U.S. Department of Energy under contract DE-AC52-06NA25396. References 1. Lucas LJ, Owhadi H, Ortiz M (2008) Rigorous verification, validation, uncertainty quantification and certification through concentration-ofmeasure inequalities. J Comput Methods Appl Mech Eng 197:4591–4609 2. Hemez FM, Atamturktur SH, Unal C (2009) Prediction with quantified uncertainty of temperature and rate dependent material behavior. In: 11th AIAA non-deterministic approaches conference, Palm Springs, May 4–7, 2009 3. Maupin R, Hylok J, Rutherford A, Anderson M (2005) Validation of a threaded assembly, part I: overview. In: 6th European conference on structural dynamics, Paris, Sept 5–7, 2005 4. Hylok J, Rutherford A, Maupin R, Anderson M, Groethe M (2005) Validation of a threaded assembly, part II: experiments. In: 6th European conference on structural dynamics, Paris, Sept 5–7, 2005 5. Rutherford A, Maupin R, Hylok J, Anderson M (2005) Validation of a threaded assembly, part III: validation. In: 6th European conference on structural dynamics, Paris, Sept 5–7, 2005 6. Owhadi H, Sullivan TJ, McKerns M, Ortiz M, Scovel C (2010) Optimal uncertainty quantification. Technical report 2010-03, Applied and Computational Mathematics, California Institute of Technology, Pasadena 7. McDiarmid C (1989) On the method of bounded differences. In: Siemons J (ed) Surveys in combinatorics. Volume 141 of the london mathematical society lecture. Cambridge University Press, Cambridge, pp 148–188 8. The MathWorks, Inc. (2011) MATLAB R2011a product help, Natick 9. Rosenbrock HH (1960) An automatic method for finding the greatest or least value of a function. Comput J 3:175–184 10. Hemez FM, Atamturktur SH (2011) The dangers of sparse sampling for the quantification of margin and uncertainty. Reliab Eng Syst Saf 96:1220–1231 11. Hemez FM (2010) Performance of the Morris one-at-a-time sensitivity analysis to explore large-dimensional functions. Technical report LAUR-10-0069, Los Alamos National Laboratory, Los Alamos 12. Stull CJ, Hemez FM (2012) Optimal inequalities to bound a performance probability using McDiarmid concentration of measure theory. Technical report LA-UR-12-24769, Los Alamos National Laboratory, Los Alamos

RkJQdWJsaXNoZXIy MTMzNzEzMQ==