A new reliability-based data-driven approach for noisy experimental data with physical constraints

A new reliability-based data-driven approach for noisy experimental data with physical constraints Jacobo Ayensa-Jim´eneza, Mohamed H. Doweidara, Jose A. Sanz-Herrerab, Manuel Doblar´e1,∗ aMechanical Engineering Department, University of Zaragoza, Spain; Arag´on Institute of Engineering Research (I3A), University of Zaragoza, Spain; Centro de Investigaci´on Biom´edica en Red en Bioingenier´ıa, Biomateriales y Nanomedicina (CIBER-BBN), Spain. Campus R´ıo Ebro, Edificio I+D, Mariano Esquillor S. N. bSchool of Engineering (ETSI), University of Seville, Spain Abstract Data Science has burst into simulation-based engineering sciences with an impressive impulse. However, data are never uncertainty-free and a suitable approach is needed to face data measurement errors and their intrinsic randomness in problems with well-established physical constraints. As in previous works, this problem is here faced by hybridizing a standard mathematical modeling approach with a new data-driven solver accounting for the phenomenological part of the problem, with the aim of finding a solution point, satisfying some constraints, that minimizes a distance to a given data-set. However, unlike such works that are established in a deterministic framework, we use the Mahalanobis distance in order to incorporate statistical second order uncertainty of data in computations, i.e. spread and correlations. We develop the underlying stochastic theoretical framework and establish the fundamental mathematical and statistical properties. The performance of the resulting reliability-based data-driven procedure performance is evaluated in a simple but illustrative unidimensional problem as well as in a more realistic solution of a 3D structural problem with a material with intrinsically random constitutive behavior as concrete. The results show, in comparison with other data-driven solvers, better ∗Corresponding author Email address: [email protected] (Manuel Doblar´e) Preprint submitted to Journal of L A T EX Templates July 21, 2017 This is an Accepted Manuscript of an article published by Elsevier in Computer Methods in Applied Mechanics and Engineering, Volume 328, on January 2018, available at: https://doi.org/10.1016/j.cma.2017.08.027 Cop. 2017 Elsevier. En idUS: Licencia Creative Commons CC BY-NC-ND convergence, higher accuracy, clearer interpretation, and major flexibility besides the relevance of allowing uncertainty management, with low computational demand. Keywords: Data-Driven, Reliability, Mahalanobis distance 1. Introduction Nowadays, Data Science and disciplines such as Big Data or Data Analytics [1] are essential in our everyday life. Photos and videos handling, control of patients data, consumer preferences data, census information and police incident reports are just some examples of the daily huge data treatment.5 These methodologies permit the extraction of patterns and/or relevant information from available unstructured data [2]. Since the main ideas and concepts were introduced at the beginning of the century [3, 4], an extensive literature may be found on this broad area [5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. For example, in Machine Learning approaches [3, 15], the idea is improving10 continuously the accuracy of predictions by means of new available data, adding new explicit knowledge from the actual response to previous predictions. A particular subdiscipline of Machine Learning is Manifold Learning [16, 17] in which the particular aim is getting newer and richer hidden knowledge related to the underlying structure or, in mathematical terms, the dimensionality and local15 bases of the relevant working space. There are many methods for building the underlying manifold from data, ranging from pure interpolation to pure regression, including all Manifold Learning techniques (kernel Principal Component Analysis (kPCA) [18], Self Organizing Map (SOM) [19], Locally Linear Embedding (LLE) [20], Isomap [21], Laplacian Eigenmap [22], t-distributed Stochastic20 Neighbor Embedding (t-SNE) [23] among others [8]). In the same direction, since Rosenblatt developed the perceptron [24], artificial neural networks is another field where new concepts as Deep Learning and dynamic networks are in continuous development and are able to identify more abstract features and solve more complex problems [25, 26, 27].25 2 Despite the wide application of Data Science in areas such as marketing and e-commerce [28], social sciences [29], or healthcare [30], there are other fields where very little has been done. An example are the disciplines where physical models and the corresponding mathematical and numerical simulation tools are well established like Computational Physics, Computational Chemistry30 or Computational Engineering (Simulation-Based Engineering and Sciences - SBES-). A straightforward application of these techniques is Dynamic DataDriven Applications Systems (DDDAS) [31], in which the idea is providing both predictive and learning capabilities to the control system of data acquired from a sufficient set of sensors. This paradigm was settled down by Kalman [32] in35 the sixties with his groundbreaking filter and is still nowadays a hot topic of research opening up a huge range of possibilities [33]. Some work has been done for dynamical systems [34] and for parameterized PDEs systems [35]. SBES may incorporate, in addition to data, some a priori characteristic physical knowledge of the analyzed system. At this point, it is crucial to distin-40 guish between two kinds of knowledge. On one hand, physical general principles, such as conservation and thermodynamic laws that are universally accepted as able to describe the underlying universe structure. On the other hand, we find phenomenological models, such as macroscopic material constitutive relations. The latter is an intelligent simplification of the real interactions at molecular45 level extracted from available experimental data. From above, it is clear that Data Analytics techniques would be very useful in SBES to extrapolate the phenomenological model, but now constrained with the mathematical expression of first principles. This approach, of increasing importance, is known as Data-Driven Simulation-Based Engineering and Sciences50 (DDSBES). In this mixed approach, the absence of physical constraints implies recovering the Data Science and Machine Learning framework, while total a priori parameterization of experimental data recovers classical SBES. Actually, all linear and non-linear phenomenological models are formulated in terms of parametric mathematical equations, where the variables of interest are forced55 to remain within a given pre-established manifold. DDSBES may be considered 3 then as an a posteriori manifold constructor that may be context-dependent. In other words: let the data tell us which physical variables persist without forcing them a priori, except for universal physical laws. Recently, F. Chinesta and coworkers defined a strategy for Data-Driven (DD)60 Computational Mechanics [36], combining Manifold Learning techniques and a (possibly optimized) directional search strategy inspired in the LaTin method [37]. In that work, it is highlighted that manifold construction step is not compulsory, but could give some insight about underlying physical structure and could result in less computationally demanding solutions. M. Ortiz and65 his group presented a material model-free method based on the minimization of the distance between the searched solution and a set of experimental data, using a proper energy norm, while remaining in the equilibrium manifold, or equivalently, a well-posed penalty approach [38]. Other DD hybrid approaches use Gaussian processes (GP) in a given data-set for dynamical feedback of the70 parametrical model [39], or fusion prognosis [40, 41], to combine DD and physical models, but these methodologies remain encapsulated in the underlying specific physical phenomenological model. None of these works take into consideration the inherent inaccuracy of the data. Conversely, empirical data are considered as error-free for both directional75 search and penalty approaches [36, 38]. Only for the latter, some mathematical convergence results are derived for zero uncertainty approaches, which in practice is never the case. In this paper, a new family of methods, called reliability-based data-driven solvers (RBDD), based on a metric accounting for uncertainty are defined and80 some new mathematical results are derived. It is highlighted how DD solver methodology naturally allows incorporating reliability along the statement of the modeling. The penalty approach suggested in [38] is chosen as starting point, but now, a metric taking into account the uncertainty, the Mahalanobis distance, is employed, in order to deal with spread and correlations of the data.85 With this, data-driven simulations become sensitive to measurements precision and incorporate uncertainty considerations. 4 Through this paper, we will discuss the main problems of the classical and the new simulation-based techniques when dealing with noisy data, highlighting the limitations of each methodology. An easy but illustrative one-dimensional prob-90 lem is used to compare results and to show improvements using this methodology. We also present another more realistic example with real concrete test data, emphasizing the implications of nonexistence of an explicit set of well-defined hypotheses and the corresponding material model. 2. Data-driven solvers95 Following [38] and [36], DD solvers may be seen as iterative solvers searching for the intersection of a (data based) empirical manifold and a physical manifold. The first one is in many practical applications experimentally based and has, therefore, a discrete nature. The second is usually established in terms of sound laws particular to the problem in hands, but otherwise derived from first100 principles universally accepted as the basis of Physics. For the sake of simplicity, we may consider the elastic three-dimensional problem. In that case, the physical manifold is the set of states that verify global and local equilibrium (i.e. conservation of linear and angular momenta), that in the static case (negligible inertial effects) is written in differential form as:105 ∇·σ=0(1) with σthe stress tensor. Equation (1) is usually approximated and solved in a discrete form using numerical methods like Finite Elements (FEM). In that case, after a convenient discretization we can state:110 By =0(2) where yis a finite dimensional vector containing the full stress tensor field information related to a given discretization (for FEM, this vector contains the 5 components of the stress tensor for all the integration points) and Bis a matrix encoding the geometry and connectivity of the domain. The empirical manifold is defined via a set E={(xj;yj)}j=1,··· ,m of data115 points, resulting from experimental measurements (and therefore not uncertainty free) as it will be illustrated subsequently. The set Emay be seen as a representation of the underlying material behavior in the following asymptotic sense: (i) if Eapproximates a mathematical manifold and (ii) uncertainty of each point approximates to zero. Some basic mathematical results related to120 these considerations may be found in [38]. When solving the problem, there are two main approaches: 1. The first one is based on identifying, at least locally, manifolds from data. Here, regression techniques (based on least squares or other optimization approach) [36] or interpolation techniques [42, 43] are generally used. For125 high dimensional spaces, regression algorithms are expensive and therefore a previous step including dimensionality reduction, i.e. Manifold Learning, is generally compulsory [36]. Once the underlying manifold is built, locally tangent spaces may be computed and tangent-based iterative solvers such as Newton-Raphson (NR), quasi-Newton or arc-length strategies may be130 used. We call this the linearization approach. 2. The second one is searching for the solution point directly from data, following [38], or equivalently the third approach presented in [36]. In that case, it is necessary to define a distance, i.e. a metric, in order to select the nearest data point to the physical manifold. In [36], the135 euclidean norm is selected, despite it is dimensionally inconsistent, while in [38] a more physically-meaningful norm (energy norm) is selected. Both norms, however, do not consider uncertainty in data. The problem is then formulated as a constrained minimization problem, and solved iteratively. We call this the pure DD approach.140 6 2.1. Problem formulation We present here the general framework for DDSBES problems. With this aim, we postulate that a model-free engineering problem may be defined in terms of state variables (X, Y ) that are related through a latent and unknown relationship F(X, Y ) = 0. For most computational frameworks, state variables145 are presented in a discrete manner such as (X, Y ) = (xi,yi)i=1,··· ,N where xi, and yiare vectors whose dimension nis the size of the state vector and N=N×nis the number of scalar state variables of the problem. Returning to the elastic problem, xiis the vector containing all strain components (εkl) at the point iand ythe vector containing all stress components (σkl) at the point150 i. It is now necessary to define a distance (a metric) in the state space for xi and yi. That is, for example for xi, to define a symmetric and positive-definite matrix Mxand: ||xi||2 x,i =1 2xT iMx,ixi(3) Therefore: d2 x,i(xi,x0 i) = ||xi−x0 i||2 x,i =1 2(xi−x0 i)TMx,i(xi−x0 i) (4) As we are considering engineering problems, we have physical constraints.155 For the sake of simplicity, but without any conceptual limitation, we shall consider linear constraints only, so they can be written as: Ax =a Cy =c(5) At each point i, we have a trial set Eithat may be thought as the result of experimental tests. We then define a local penalty function for each point ias: Fi(xi,yi) = min (x0,y0)∈Ei {dx,i(xi,x0) + dy,i(yi,y0)}(6) 7 It is obvious that the penalty function vanishes for each point in Ei,Fi|Ei= 0.160 Finally, a global penalty function is defined, F(x,y|E) = PN i=1 Fi(xi,yi) where E=QN i=1 Ei,x= (xi)i=1,··· ,N and y= (yi)i=1,··· ,N . Here we have, for each (x,y)∈ E,F(x,y) = 0, F|E= 0 and we have a global norm ||(x,y)||2= PN i=1 ||(xi,yi)||2. Therefore, a DDSBES problem is defined by the constrained optimization165 problem: min (x,y)F(x,y|E) subject to Ax =a Cy =c (7) For the elastic problem, problem (7) takes the form: min (ε1,··· ,εN,σ1,··· ,σN)F(ε1,· · · ,εN,σ1,· · · ,σN|E) subject to C[σ1,· · · ,σN]T=c (8) where Cis a matrix encoding connectivity and geometry of the problem, depending on the particular discretization. This formulation is similar to the one proposed in [38], except for the fact170 that it is formulated in a slightly more general context, including a generalized distance (and therefore a more flexible way of measuring how far is a point from the data-set). It is important to note that state variables xi,yimay be in a lower dimensional space obtained after dimensionality reduction. For example, for the175 elastic problem, εand σlive in a space of 6 dimensions, so (ε,σ) has 12 dimensions. These dimensions may be reduced if additional simplifications are imposed a priori onto the material behavior. For example, if a homogeneous isotropic linear material is considered, this dimension is actually two, because 8 of the Hooke law states σ=λTr(ε)+2µεwhere λand µare Lame’s constants,180 related to phenomenological parameters Eand ν. Proposition 1. The problem defined by (7) has a unique solution if rang(A) = rang(C) = rwhere ris the number of restrictions. Proof:. Let Mx=Ln i=1 Mx,i yMy=Ln i=1 Mx,i. Let (x∗,y∗) a pair of state variables verifying Fi(xi,yi) = d2 x,i(xi,x∗ i)+d2 y,i(xi,y∗ i) which exists because of185 Eifiniteness. With these definitions we can write F(x,y) = 1 2(x−x∗)TMx(x−x∗) + 1 2(y−y∗)TMy(y−y∗) (9) We define the lagrangian function L(x,y,λ,µ) = F(x,y)−λT(Ax −a)− µT(By −b). Then: ∂L ∂x=Mx(x−x∗)−ATλ ∂L ∂y=My(y−y∗)−BTµ ∂L ∂λ=Ax −a ∂L ∂µ=By −b(10) Using Lagrange multipliers theorem: ∂L ∂x=0,∂L ∂y=0,∂L ∂λ=0y∂L ∂µ=0, therefore:190 Mxx−ATλ=Mxx∗ Myy−BTµ=Myy∗ Ax =a By =b(11) 9 Figure 1: Rod under uniaxial load. It is important to note in the next figures, corresponding to the numerical experiments, that the grey dashed line, drawn as a function graph σ=f(), represents the actual material behavior obtained in the experiments, whose expression is not known a priori. When using linearization approaches, we need to310 define, at least locally, a smooth manifold in order to work with tangent spaces. This can be done by using the multiple Manifold Learning techniques presented at the Introduction. As a rule of thumb, the more accurate and structured empirical data set, the better interpolation-based techniques perform. On the other hand, regression techniques are preferred when dealing with noisy and315 unstructured data but low dimensional and regular underlying manifolds are desirable. Besides, Ecan be generated either allowing control in one of the variables (laboratory controlled tests) or control is impossible (for example sensors in dynamic DD systems). The later case is the most general and challenging. We320 are going to test the convergence for these two cases using the typical NewtonRaphson solver. We fixed a maximum number of iterations to 104which is huge taking into consideration that, usually, this kind of solvers achieve convergence in a few iterations. Four analyses are considered, varying the number of data sample points, m,325 and the error measure related to uncertainty, s. Data generation is as follows: for each εj∼ U(0; εmax), j= 1,· · · , m and, as before, σj∼ N (µj, s), with 16 µj=f(εj) and s=ασmax, where σmax is the maximum stress. Figures 2a, 3a, 4a, 4c, 5a and 5c show the considered empirical set, the equilibrium manifold and the constitutive manifold built for some fitting tech-330 niques (linear interpolation, natural spline interpolation and 5 degree polynomial regression). The vertical dashed line shows initial points considered for the Newton-Raphson solver. Figures 2b, 3b, 4b, 4d, 5b and 5d show the empirical set, the equilibrium manifold and final point for each solver. Both reliabilitybased data-driven (RBDD) and DD solvers converge to the same point. Con-335 vergence is not achieved by the Newton-Raphson solver based on regression fit because of the untrue local convexity of the built manifold, which is inherent to parametric regression. In the case of natural spline regression, almost linear behavior in the hardening part of the curve causes bad convergence. In Figure 2, solvers are based on an empirical set of m= 100 pairs (ε;σ) and a low340 homogeneous uncertainty is considered (α= 0.001). 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Linear interpolation Natural spline interpolation 5 degree polynomial regression Equilibrium manifold Initial point (a) Built manifold for linearization techniques. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Equilibrium manifold Linear interpolation + NR Natural spline interpolation + NR Regression + NR DD & RBDD solvers (b) Solution point for different solvers. Figure 2: Performance of different solvers for m= 100 and α= 0.001. In Figure 3, we use an empirical set of m= 100 pairs (ε;σ) but higher homogeneous uncertainty is considered (α= 0.05). Even if polynomial regression is not sensitive to noise, convergence is again not achieved because of the untrue local convexity of the built manifold. Besides, due to noise, natural splines suf-345 fer spurious oscillations provoking bad convergence. This can be avoided using 17 linear interpolation, but in this case, non-smoothness of the broken line is incompatible with a tangent-based solver, which in turns results in non-convergence. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Linear interpolation Natural spline interpolation 5 degree polynomial regression Equilibrium manifold Initial point (a) Built manifold for linearization techniques. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Equilibrium manifold Linear interpolation + NR Natural spline interpolation + NR Regression + NR DD & RBDD solvers (b) Solution point for different solvers. Figure 3: Performance of different solvers for m= 100 and α= 0.05. We next analyze the solver behavior for reduced sample sizes. Using an empirical set of m= 20 pairs (ε, σ). First, we consider accurate data (α= 0.001).350 For soundness considerations, we analyze the case with F= 100 kN and F= 200 kN. Empirical sets are different but are associated to the same m and α. Results are shown in Figure 4. The regression based solver is not convergent for any method. Obviously, due to the lack of data, the DD solver has, in that case, less accuracy than linearization approaches based on interpola-355 tion techniques. However, spline interpolation may also have bad convergence, depending on the empirical set mesh and the equilibrium manifold. Linear interpolation would give accurate results only for quasi-linear behavior and/or fine constitutive manifold meshes. Convergence problems increase dramatically when considering greater noise, as seen in Figure 5, where only DD solvers360 converge to an accurate enough solution. For homogeneous uncertainty, the DD solver and the RBDD solver give the same result, as pointed out before. Table 1 shows the squared distance results for DD solvers. RBDD is more informative in the following sense: for m= 100, the distance to the empirical set increases when passing from α= 0.001 to365 18 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Linear interpolation Natural spline interpolation 5 degree polynomial regression Equilibrium manifold Initial point (a) Built manifold for linearization techniques, F= 100 kN. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Equilibrium manifold Linear interpolation + NR Natural spline interpolation + NR Regression + NR DD & RBDD solvers (b) Solution point for different solvers, F= 100 kN. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Linear interpolation Natural spline interpolation 5 degree polynomial regression Equilibrium manifold Initial point (c) Built manifold for linearization techniques, F= 200 kN. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Equilibrium manifold Linear interpolation + NR Natural spline interpolation + NR Regression + NR DD & RBDD solvers (d) Solution point for different solvers, F= 200 kN. Figure 4: Performance of different solvers for m= 20 and α= 0.001. 19 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Linear interpolation Natural spline interpolation 5 degree polynomial regression Equilibrium manifold Initial point (a) Built manifold for linearization techniques, F= 100 kN. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Equilibrium manifold Linear interpolation + NR Natural spline interpolation + NR Regression + NR DD & RBDD solvers (b) Solution point for different solvers, F= 100 kN. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Linear interpolation Natural spline interpolation 5 degree polynomial regression Equilibrium manifold Initial point (c) Built manifold for linearization techniques, F= 200 kN. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data True behaviour Equilibrium manifold Linear interpolation + NR Natural spline interpolation + NR Regression + NR DD & RBDD solvers (d) Solution point for different solvers, F= 200 kN. Figure 5: Performance of different solvers for m= 20 and α= 0.05. 20 α= 0.05. This is due to pure hazard; each realization will give us a different distance depending only on the empirical set sample. RBDD solver does not have this problem because it is uncertainty dependent and can detect when uncertainty is of the order of the optimal distance. Only when α→0 and the empirical set is almost a subset of the constitutive manifold, this distance370 can be used as a good uncertainty-free indicator. Otherwise, the locus of the underlying manifold is unknown and there is no way to interpret DD optimal distance in a coherent manner. m100 20 20 F [kN] 100 100 200 α0.001 0.05 0.001 0.05 0.001 0.05 DD solver 275 6040 75251 112431 6692 1602 RBDD solver 183.43 1.61 50167.62 29.98 29.59 0.096 Table 1: Squared distance results for DD solvers. To analyze the statistical properties of the squared distance, we consider α= 0.1, m= 100 and F= 230 kN to generate two possible empirical sets with375 the same statistical properties. Let us suppose that the expected value of the empirical set is known. This should approximate the true underlying manifold but actually it may be estimated from experimental samples. Figure 6 shows the considered data points, expected values µand an error band, defined by µ±s, where sis the standard deviation. Table 2 shows statistical properties of380 the RBDD solver for both cases, assuming normality. It is important to note that here we consider εas an uncertainty free variable, and therefore chi-squared distribution of D2 E=σ−µ s2has n= 1 degree of freedom. The squared distance computed from data is almost zero in both cases, but a deeper knowledge about empirical set statistics (µ=µ(ε), s=s(ε)) highlights the distance to the true385 manifold. Anyway, for the case analyzed, 96 simulations have had to be carried out to obtain such result (Case 2). RBDD solver is not only a more suitable and more informative solver. It 21 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data DD & RBDD solution Equilibrium manifold 7 7 ' s (a) Solution near to the underlying manifold (Case 1). 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 2 4 6 8 10 12 14 16 18 <[MPa] Data DD & RBDD solution Equilibrium manifold 7 7 ' s (b) Solution far from the underlying manifold (Case 2). Figure 6: Performance of DD and RBDD solvers for m= 100 and α= 0.1 considering two different samples. Case 1 2 Squared Optimal distance 0.0056 0.0000 Expected value 1.001 6.725 Variance 2.003 24.899 95%-Confident bound 3.84 16.30 Table 2: Statistical characteristics of the two solution points. 22 can, for non-homogeneous uncertainty, result in a proper convergence in the following sense. Figure 7 shows solution points for the DD and RBDD solvers390 for a material with different (ε, σ) constitutive relationship. Note that the uncertainty associated with the actual material behavior is not homogeneous: in the elastic zone, where the material is very well characterized, uncertainty is low, but it increases when strains are higher. RBDD solver is sensitive to this variation, while DD solver is not. For complete information, Figure 7 should395 be complemented by the statistical properties summarized in Table 3. Thus, in Figure 7a, we can see that the DD solver converges to a very unlikely point while RBDD converges to a more likely one. This is due to the smaller ratio between the geometric distance from the solution point to the empirical data set point and the bandwidth, in case of the RBDD. However, even though the400 squared distance is small, the expected value indicates that the RBDD solver has converged to a point not very close to the mean manifold. In Figure 7b, a different convergence point is also observed. Now, the RBDD solver has an undesirable behavior because of the lack of data in the linear and certain region. This is detected by means of the expected value and variance, as well as the405 95% upper bound. RBDD solvers may be therefore used for sampling strategy considerations. In any case, this is a very unreasonable case, because often, more sample points are associated with less uncertainty. In Figure 7c we can see that both solvers converge to the same values and the statistical properties are similar to those of the first case, indicating a reliable convergence. Finally, the410 fourth case is similar to the second one, but with higher uncertainty, which in turn reduces the expected value of the squared distance, although the variance remains relatively high. Note that knowledge of the upper confident bound of the squared distance, D2 E, could be interesting for defining a quantitative criterion for convergence.415 3.2. Scale data reduction An interesting application of the here introduced RBDD solver in the domain of Computational Mechanics appears when dealing with several scales. 23 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 5 10 15 <[MPa] Experimental data RBDD solution DD solution Equilibrium manifold 7 7 ' s (a) Case 1. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 5 10 15 <[MPa] Experimental data RBDD solution DD solution Equilibrium manifold 7 7 ' s (b) Case 2. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 5 10 15 <[MPa] Experimental data RBDD solution DD solution Equilibrium manifold 7 7 ' s (c) Case 3. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0 5 10 15 <[MPa] Experimental data RBDD solution DD solution Equilibrium manifold 7 7 ' s (d) Case 4. Figure 7: Performance of DD and RBDD solvers for heterogeneous uncertainty. Case 1 2 3 4 DD Squared distance 24.72 ·10254.14 ·1026.98 ·10254.14 ·102 RBDD Squared distance 0.08 1.12 1.44 0.01 Expected value 3.70 10.98 4.20 4.94 Variance 12.82 41.91 14.18 17.76 95%-Confident bound 10.82 23.07 11.79 13.18 Table 3: Statistical properties of both solvers for each of the presented cases. 24 One of the main strategies used when coupling two scales (multiscale approach) is selecting a representative volume element (RVE) and establishing a sound420 transition procedure between the microscale properties and the macroscale response [45, 46]. This strategy has allowed setting up implicitly material constitutive relationships that were not known explicitly at the macroscale [47, 48]. Recent works foreground the crucial point of scales decoupling in the averaging process and the need of uncertainty quantification when building the restriction425 operator [49]. As answer, many works have incorporated microscale randomness in the multiscale procedure, either using Montecarlo Method (MCM) sampling [50, 51, 52] or Perturbation Method [53]. However, these considerations are still used for model validation and uncertainty has not been incorporated routinely in macroscale computations, except through expensive MCM sampling.430 A different alternative has been proposed by using stochastic partial differential equations (SPDEs) [54]. We can apply the presented RBDD solver for uncertainty propagation from the microscale to the macroscale allowing to incorporate it in macroscale computations. Let us assume that we have at the microscopic scale a (discrete)435 coupled field (xi,yi), i = 1,· · · K. Therefore, classical RVE techniques allow us to define a macroscopic reference value (X,Y) where X=xiand Y=yi. For elastic problems, this could be strain-stress pairs (εi,σi). It is possible to compute the variance-covariance matrix Σof the sample ui= (xi,yi). Geometrically, this means to define a 2n-dimensional ellipsoid in the state space (X, Y )440 associated to each single macroscopic point, where nis the space dimension of state variable xor y. Figure 8 shows an ellipsoid in a two-dimensional plane for two possible microstructural fields. Note that accounting only for average values, as done in classical RVE techniques, gives the same result in both of them. We have presented here an approach from the point of view of dimensionality re-445 duction (we use first and second order statistics instead of the whole microscopic field). Besides, this could be exploited in multiscale computational mechanics if a resourceful method allows second order statistical characterization of the microscale without the whole microscale fields computation. 25 Figure 13: Experimental data, confidence band at level µ±sand numerical solution ad different points and steps. 4. Discussion and conclusions In this work, a new RBDD solver has been formulated for DDSBES problems, allowing uncertainty considerations in the input data that are, therefore, not considered as uncertainty-free, but of random nature. The DDSBES problem is here defined as a constrained stochastic optimization problem. Constraints535 encode all relevant physical information of the system, such as fundamental conservation or physical laws. Calculations are carried out directly from data, avoiding any modeling error via state or constitutive equation assumptions. Optimality is sought in terms of a penalty function showing the distance between a candidate solution point and input data set.540 It has been shown that employing a proper uncertainty dependent distance, the Mahalanobis distance, results in good statistical properties as well as an easily interpretable optimal distances. Indeed, this optimal distance is computed in terms of the data sample and data uncertainty, which allows assessing when the solution point is accurate enough up to data precision. Moreover, this545 distance offers the possibility of considering heterogeneous uncertainty, leading to most likely solution points, instead of getting deterministic solution points which may be very sampling-dependent. 32 Excluding very simple problems with uncertainty-free linear behavior and small sample sizes (Figure 4 and Figure 7), where conventional solvers could be550 used, the method here proposed has shown better convergence, higher precision, clearer interpretation, major flexibility and more soundness. Besides, the statistical interpretation, depending on sampling statistics, allows decision-oriented statistical inference. RBDD solvers appear to be a very suitable tool for facing at least three555 important problems: 1. Dynamic DD systems, where predicting features and learning capabilities are combined. The presented solver could start predictions from scratch, where sample sets are small and the underlying true manifold is unknown. Further knowledge of the analyzed system due to the increase of the sam-560 ple size will feed the RBDD solver, thus allowing faster updates of solution points and statistical inference. This will, additionally, enhance the possibility to define different sampling strategies, data coverage and improve solver performance. Moreover, it is sensitive to measurement errors, that depends on equipment and human precision. In this sense, RBDD solvers565 conform a robust framework that provides coherent results within the experimental context. 2. Scale dimensionality reduction problems. From purely theoretical (sound physical reasons) and/or practical (speed-up calculations) point of views, it could be interesting to define the transition from a small scale to a higher570 one, defining a hierarchical procedure. The presented RBDD is an ideal tool for uncertainty propagation from one scale to another. Moreover, this may be helpful in case of a dimensionality reduction strategy anchored to Big Data frameworks. For instance, if we work with an n-dimensional field and two scales with two mesh sizes Nand M, the whole problem will575 have N·M·ndegrees of freedom. Averaging techniques could reduce the problem to a M·ndegrees of freedom problem, but all the variability of the lower scale is lost. With the RBDD approach, lower scales variability 33 is conserved, at low computational cost, resulting in M·n(n+3) 2∼θ(Mn2) degrees of freedom. If nM, savings are evident.580 3. Model-free engineering based on empirical measurements. The new presented solver offers the possibility of carrying out simulation directly from data, without explicit model assumptions. However, this rises an important limitation not only of this methodology, but of Data Analytics in general. This corresponds to the need of contextualize data in order to585 be sure that they can be extrapolated to a possibly different context corresponding to the particular application. The possibility (or not) of this extrapolation is analyzed explicitly in the standard model-driven approach when making explicit the assumptions that drive to such particular model. The need of matching the DD methodology with existing (simplified) exper-590 imental setups today available to capture the mechanical behavior of materials implies the need of making some explicit assumptions or, at least, to think about the context in which the data have been obtained and the one of the application in hand to decide if they can be extrapolated, and if there is additional data required to fulfill the problems demands. Moreover, this method relies on the595 hypothesis that we have complete information for each point of the data-set, that is, for each point, all the coordinates are known (for example, in the problem arising from computational mechanics, all the components of both tensors ε and σare known for each data point). When this is not the case, an appropriate filling data strategy should be considered.600 Regardless of this, RBDD solvers present a meeting point between theoretical sciences, through epistemologic constraints, and experimental sciences, through uncertain real world data. The elegance of the mathematical formulation enables many analysis and theoretical considerations for the whole spectrum of Continuum Physics. The ease of combining the presented concepts with all605 trendy Data Science and Deep Learning tools opens up huge possibilities for facing the most challenging problems because it offers a huge range of possibilities in dynamic DD applications, dimensionality reduction, decision-support 34 systems and any kind of problem in which uncertainty plays a major role. Acknowledgments610 The authors gratefully acknowledge the financial support from the Spanish Ministry of Economy and Competitiveness (MINECO MAT2016-76039-C4-4R, AEI/FEDER, UE), the Government of Aragon (DGA) and the Biomedical Research Networking Center in Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN). CIBER-BBN is financed by the Instituto de Salud Carlos III615 with assistance from the European Regional Development Fund. 35 Appendix A. Mathematical proofs Definition:. Let X= (X1,· · · , Xn)Ta random vector and Ma symmetric positive-definite matrix. The stochastic quadratic form (SQF) QM(X)is the random variable defined620 as: QM(X) = XTMX (Appendix A.1) Lemma Appendix A.1. Let QA(X)and QB(X)two SQF, and let {ek}k=1,··· ,n the standard basis in Rn, then: E[QA(X)] = AijΩij(X) (Appendix A.2) E[QA(X)X] = AijΛijk(X)ek(Appendix A.3) E[QA(X)QB(X)] = AijBklΥijkl(X) (Appendix A.4) Cov(QA(X), QB(X)) = AijBkl (Υijkl(X)−Ωij(X)Ωkl(X)) (Appendix A.5) where:625 Ωij(X) = E[XiXj] (Appendix A.6) Λijk(X) = E[XiXjXk] (Appendix A.7) Υijkl(X) = E[XiXjXkXl] (Appendix A.8) 36 Proof:. Using index notation QA(X) = AijXiXjtherefore: E[QA(X)] = E[AijXiXj] = AijE[XiXj] = AijΩij(X) (Appendix A.9) E[QA(X)X] = E[AijXiXjXkek] = AijekE[XiXjXk] = AijΛijk(X)ek (Appendix A.10) E[QA(X)QB(X)] = E[AijXiXjBklXkXl] = AijBklE[XiXjXkXl] = AijBklΥijkl(X) (Appendix A.11) Finally, Cov(QA(X), QB(X)) =E[(AijXiXj−AijE[XiXj])(BklXkXl−BklE[XkXl])] =AijBklE[XiXjXkXl−E[XiXj]XkXl−XiXjE[XkXl] + E[XiXj]E[XkXl]] =AijBkl (E[XiXjXkXl]−E[XiXj]E[XkXl]−E[XiXj]E[XkXl] + E[XiXj]E[XkXl]) =AijBkl (Υijkl(X)−Ωij(X)Ωkl(X)) (Appendix A.12)  Proposition Appendix A.1 (Expectation of a SQF). Let QM(X)a SQF and let µ(X)the expected value of Xand Σ(X)the variance - covariance matrix630 of X. Therefore: E[QM(X)] = Tr(MΣ(X)) + µ(X)TMµ(X) (Appendix A.13) Proof:. Following Lemma Appendix A.1, E[QM(X)] = MijΩij. Furthermore, Ωij(X) = Σij(X) + µi(X)µj(X) then: 37 E[QM(X)] = Mij(Σij(X) + µi(X)µj(X)) =MijΣij(X) + Mijµi(X)µj(X) = Tr(MΣ(X)) + µ(X)TMµ(X) (Appendix A.14) In the last equality we have used Tr(AB) = A:B. 635 Proposition Appendix A.2 (Variance and covarianze of SQF). Let QA(X) and QB(X)two SQF and let µ(X)the expected value of X,Σ(X)the variance - covariance matrix of Xand Υ(X)the fourth order moment tensor of X. Then: Cov(QA(X), QB(X)) =A:Υ(X) : B−Tr(AΣ(X)) + µ(X)TAµ(X)Tr(BΣ(X)) + µ(X)TBµ(X) (Appendix A.15) In particular, if A=B: Var(QA(X)) = A:Υ(X) : A−(Tr(AΣ(X)) + µ(X)TAµ(X))2 (Appendix A.16) Proof:. Following Lemma Appendix A.1, Cov(QA(X), QA(X)) = AijBklΥijkl(X)+640 AijBklΩij(X)Ωkl(X). However, Ωij (X)=Σij(X) + µi(X)µj(X) therefore: Cov(QA(X), QA(X)) =AijBklΥijkl(X) + AijBkl(Σij(X) + µi(X)µj(X))(Σkl(X) + µk(X)µl(X)) =AijBklΥijkl(X)+(AijΣij(X) + Aijµi(X)µj(X))(BklΣkl(X) + Bklµk(X)µl(X)) =A:Υ(X) : B−(A:Σ(X) + µ(X)TAµ(X))(B:Σ(X) + µ(X)TBµ(X)) (Appendix A.17) The final result is obtained noting that Tr(AB) = A:B. 38  Let’s now assume normality. The following result may be found in [55]: Remark Appendix A.1 (Fourth order moments of centered multivariate normal distribution).645 Let Z∼ N(0,Σ)n-dimensional multivariate normally distributed random vector whose expected value is zero and variance - covariance matrix is Σ. Then: µi(Z) = 0 (Appendix A.18) Ωij(Z)=Σij (Appendix A.19) Λijk(Z) = 0 (Appendix A.20) Υijkl(Z)=ΣijΣkl + ΣikΣjl + ΣjkΣil (Appendix A.21) Lemma Appendix A.2. Let Z∼ N (0,Σ)an n-dimensional multivariate normally distributed random vector with expected value µ=0and variance - covariance matrix Σ. Therefore, for symmetric matrices Aand B:650 E[QA(Z)] = Tr(AΣ) (Appendix A.22) E[QA(Z)QB(Z)] = Tr(AΣ)Tr(BΣ) + 2Tr(AΣBΣ) (Appendix A.23) E[QA(Z)Z] = 0(Appendix A.24) Proof:. The first equation is obtained directly from linearity of the expected value operator and the fact that µ=0. For the second, note that, following Lemma Appendix A.1 E[QA(Z)QB(Z)] = AijBklΥijkl(Z), but, by virtue of Observation Appendix A.1, Υijkl(Z)=ΣijΣkl + ΣikΣjl + ΣjkΣil and then, using Σsymmetry:655 39 E[QA(Z)QB(Z)] = AijΣijBklΣkl +AijΣkiBklΣlj +AijΣjkBklΣli = Tr(AΣ)Tr(BΣ) + 2Tr(AΣBΣ) (Appendix A.25) Finally, for the third one, following again Lemma Appendix A.1, E[QA(Z)Z] = AijΛijkek, and then, using Observation Appendix A.1 we obtain that E[QA(Z)Z] = 0.  We can then prove the following result:660 Proposition Appendix A.3 (Variance and covariance of two SQF under normality). Let Z∼ N(µ,Σ)an n-dimensional multivariate normally distributed random vector with expected value µ=0and variance - covariance matrix Σ. Then, if Aand Bare symmetric: Cov(QA(X), QB(X)) = 2Tr(AΣBΣ)+4µTAΣBµ(Appendix A.26) In particular:665 Var(QA(X)) = 2Tr(AΣAΣ)+4µTAΣAµ(Appendix A.27) Proof:. We use the expression Cov(X, Y ) = E[XY ]−E[X]E[Y] (Appendix A.28) With X=QA(X) and Y=QB(X). The first term of the right hand side may be developed in terms of Z=X−µusing symmetry of matrices Aand B as: 40 E[QA(X)QB(X)] =E[(ZTAZ +ZTAµ+µTAZ +µTAµ)(ZTBZ +ZTBµ+µTBZ +µTBµ)] =E[ZTAZZTBZ]+2µTE[ZTAZZ]+2µTE[ZTBZZ]+4µTAE[ZZT]Bµ +E[ZTAZ]µTAµ+ +E[ZTBZ]µTBµ+ 2µTAE[Z]+2µTBE[Z] + µTAµµTBµ =E[QA(Z)QB(Z)] + 2µTE[QA(Z)Z]+2µTE[QB(Z)Z]+4µTAΣBµ +E[QA(Z)]µTAµ+E[QB(Z)]µTBµ+µTAµµTBµ = Tr(AΣ)Tr(BΣ) + 2Tr(AΣBΣ) + Tr(AΣ)µTAµ+ Tr(BΣ)µTBµ + 4µTAΣBµ+µTAµµTBµ (Appendix A.29) In last equality we have used Lemma Appendix A.2.670 The second term of the right hand side, is obtained analogously: E[QA(X)]E[QB(X)] =E[ZTAZ +ZTAµ+µTAZ +µTAµ]E[ZTBZ +ZTBµ+µTBZ +µTBµ] = (E[ZTAZ]+2µTAE[Z] + µTAµ)(E[ZTBZ]+2µTBE[Z] + µTBµ) = (E[QA(Z] + µTAµ)(E[QB(Z] + µTBµ) = (Tr(AΣ) + µTAµ)(Tr(BΣ) + µTBµ) (Appendix A.30) Again, in the last equality, Lemma Appendix A.2 was used. Subtracting Equation (Appendix A.30) to Equation (Appendix A.29), we obtain: E[QA(X)QB(X)] −E[QA(X)]E[QB(X)] = 2Tr(AΣBΣ)+4µTAΣBµ (Appendix A.31) 675 41 Highlights •A hybrid mathematical modeling and data-driven approach is proposed •Data uncertainty managed by a proper metric, the Mahalanobis distance •Good convergence, accuracy and flexibility with low computational demand were found Accepted Manuscript A new reliability-based data-driven approach for noisy experimental data with physical constraints Jacobo Ayensa-Jim´ enez, Mohamed H. Doweidar, Jose A. Sanz-Herrera, Manuel Doblar´ e PII: S0045-7825(17)30425-5 DOI: http://dx.doi.org/10.1016/j.cma.2017.08.027 Reference: CMA 11568 To appear in: Comput. Methods Appl. Mech. Engrg. Received date : 13 March 2017 Revised date : 21 July 2017 Accepted date : 20 August 2017 Please cite this article as: J. Ayensa-Jim´ enez, M.H. Doweidar, J.A. Sanz-Herrera, M. Doblar´ e, A new reliability-based data-driven approach for noisy experimental data with physical constraints, Comput. Methods Appl. Mech. Engrg. (2017), http://dx.doi.org/10.1016/j.cma.2017.08.027 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.