An unsupervised data completion method for physically-based data-driven models

Accepted Manuscript An unsupervised data completion method for physically-based data-driven models Jacobo Ayensa-Jim´ enez, Mohamed H. Doweidar, Jose A. Sanz-Herrera, Manuel Doblar´ e PII: S0045-7825(18)30488-2 DOI: https://doi.org/10.1016/j.cma.2018.09.035 Reference: CMA 12095 To appear in: Comput. Methods Appl. Mech. Engrg. Received date : 22 December 2017 Revised date : 21 September 2018 Accepted date : 24 September 2018 Please cite this article as: J. Ayensa-Jim´ enez, et al., An unsupervised data completion method for physically-based data-driven models, Comput. Methods Appl. Mech. Engrg. (2018), https://doi.org/10.1016/j.cma.2018.09.035 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. This is an Accepted Manuscript of an article published by Elsevier in Computer Methods in Applied Mechanics and Engineering, Volume 344, on February 2019, available at: https://doi.org/10.1016/j.cma.2018.09.035 Cop. 2028 Elsevier. En idUS: Licencia Creative Commons CC BY-NC-NDTu texto aquí 1 An unsupervised data completion method for physically-based data-driven models Jacobo Ayensa-Jim´eneza, Mohamed H. Doweidara, Jose A. Sanz-Herrerab, Manuel Doblar´ea,∗ 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., 50018 Zaragoza (Spain) bSchool of Engineering (ETSI), University of Seville, Spain Abstract Data-driven methods are an innovative model-free approach for engineering and sciences, still in process of maturation. The idea behind is the combination of data analytics techniques, to handle the huge amount of data derived from continuous monitoring or experimental measurements, and of the constraints imposed by universal physical laws, particular to the field in hands. A well-known problem in the former corresponds to the quality and completeness of the available data that, sometimes, are so poor that make the predictions useless. In data-driven simulation-based engineering and sciences (DDSBES), the intrinsic physical constraints may help in completing the missing data in a more precise manner, by forcing them to remain in the manifold defined by the physical laws. In this work, a suitable imputation method to complete incomplete data that preserves the data context-dependent structure is presented. This is accomplished by enforcing the set of physical constraints, specific to the problem. For this purpose, a generalization of the weighted mean concept is proposed, where the distance to the admissible points (in a physical sense) is used as a weighting function to get the optimal candidate. The method is evaluated in a classical regression problem, where it is compared with other standard methods, showing ∗Corresponding author Email address: [email protected] (Manuel Doblar´e) Preprint submitted to Journal of L A T EX Templates September 19, 2018 p2_uldd.pdf better results. Then, its application is illustrated in two data-driven problems, where no filling data procedure has been yet proposed, showing good predictive capability, provided that the data are close enough to the actual system state. Keywords: Data-driven methods, Data completion, Statistical imputation, Weighted mean, Computational Mechanics 1. Introduction Data are everywhere around us. An incredibly huge amount of sensors and transducers get measurements from the physical world. Businesses of every kind search and collect data across the globe related to consumer preferences and trends. Governments regularly collect all sorts of data from census information5 to incident reports in police departments. According to the 2016 IDC directives presented in its yearly event in San Jose (US), this deluge of data is set to rise steeply from the estimated world total amount of 4,4 zettabytes of data in 2013 to 180 zettabytes by 2025 (one zettabyte is equivalent to one trillion of gigabytes). The advent of the Internet of Things will likely make to surpass10 these figures by far [1]. More and more complicated strategies are used to extract patterns and/or relevant knowledge from this massive amount of available structured and unstructured data. In fact, the framework, in which it is easier to get predictions directly derived from available data than from tedious, complicated and some-15 times inaccurate mathematical models, is progressively changing the paradigm of predictive Physics [2, 3]. Since the main ideas and concepts were introduced at the beginning of the century, an extensive literature may be found on this broad area of Data Analytics and Artificial Intelligence [4, 5, 6]. The main trend today is the constant20 improvement of the accuracy of predictions by continuous “learning” from a non-stop input of new data, thus progressively refining the predictions by comparing the predicted and actual responses. Since the seminal idea of the perceptron [7], artificial neural networks has 2 been another pushful field where new concepts as Deep Learning and Dynamic25 Networks are in continuous development. Today, these methods allow extracting abstract features and solving very complex problems, many times not fully formalized [8, 9]. These techniques try to mimic the process of human knowledge acquisition and structuring and have become amenable after remarkable advances in sensoring; data acquisition, transfer, storage and management; enor-30 mous improvements in the performance of computers; and continuous contributions in their theoretical and algorithmic foundations. In an engineering context, a straightforward application of all these techniques is the so-called dynamic data-driven assimilation systems (DDDAS) [10], in which the idea is providing both predictive and learning capabilities to a con-35 trol system from data acquired from a set of sensors. This paradigm was settled down by Kalman [11] in the sixties with his groundbreaking filter. Nowadays, it is still a hot topic of research [12]. In the last years, a new approach to simulation-based engineering and sciences that uses the power of data-science methods has been proposed. This40 approach, of increasing importance, and known as data-driven simulation-based engineering and sciences (DDSBES), combines physical constraints and raw data. In the absence of physical constraints, the standard Data Science and Machine Learning framework is recovered, while the use of an a priori parametric model, that fits the experimental data, recovers the classical SBES. Actually,45 all linear and nonlinear phenomenological constitutive models can be formulated in terms of parametric mathematical equations, where the variables of interest are forced to remain within a given pre-established manifold. This manifold is derived from observation and experience (empiricism) eventually by means of a trial-error fitting procedure.50 One idea in this direction was started by Chinesta and coworkers [13] who defined a strategy for data-driven Computational Mechanics, combining manifold learning techniques and a (possibly optimized) directional search strategy inspired in the LaTin method [14]. Ortiz and his group [15] presented a modelfree method based on the minimization of the distance between the searched55 3 solution and a set of experimental data, using a proper energy norm. The solution was also forced to remain in the equilibrium manifold, by means of a well-posed penalty approach. This work was extended by several groups to take into account the uncertainty of the data in what are now called as reliability based data-driven solvers [16, 17].60 However, the referred works are often restricted to the frame of perfect information, that is, the data-set is complete. In other words, all state variables are assumed to be known for a given measure. This is not always the case in practical situations, being this one of the main concerns in Data Science [18, 19, 20, 21, 22, 23].65 There are many methods that have been developed to address this problem, both model-free or model-based. Among the model-free, the most fundamental are Listwise Deletion (LD) and Pairwise Deletion (PD) [22] that consist of discarding incomplete data. These methods, however, decrease statistical power [24] and introduce bias [25] if the missing process is Missing Not At Random70 (MNAR) [18, 26, 27], what is obviously the case when dealing with data obtained by experiments or measurements. Another common approach is Single Imputation (SI), based on a filling strategy for the missing data that uses values obtained from complete data, for example, Mean or Mode Imputation [19, 20]. This technique, however, reduces the variability and weakens the covariance be-75 tween variables. Another approach is to create dummy variables accounting for the missing data variables (Dummy Variable Adjustment) [19, 20]. This results in biased estimators and is not theoretically based. Finally, it is also possible to replace missing values with a predicted score from a regression equation [20, 19]. This weakens the variance and overestimates model fit and correlation80 estimates. Moreover, these methods do not take into account the local structure and geometry of the data, which is critical when the data have some underlying physics. In order to solve this problem, interpolation (linear interpolation, nearest interpolation or spline interpolation, [28]) is a common technique. In this approach, only the physics inherent to the data is learned in the imputation85 process. 4 Other model-based methods have been developed to deal with the missing data problem, such as Multiple Imputation [29] using the regression method [30], the Predictive Mean Matching Method [31] or the Markov Chain Montecarlo Method (MCMC) [25, 32], Full Information Maximum-Likelihood (FIML)90 estimation [33, 34, 35] and Expectation-maximization [36]. The problem of all these model-based methods is that they assume, one way or the other, a statistical model for the data (e.g. normality). This is usually the case in social and economical sciences [37] but is not the general case for physical problems, where variables follow some fundamental laws incompatible with normality or95 other distributional assumptions. The amount of missing data is not, however, the sole criterion to assess the quality of the available data, especially if they correspond to a problem that relies on some physical laws [38]. Our aim in this work is then to establish a framework in which both the local structure of the data and the supplementary100 physics, not explicitly included in the data structure, are used to improve the imputation procedure. In this context, the imputation method can be compatible with any DDSBES method. The presented technique is based on the mean concept and, therefore, could be interpreted as a generalization of the Mean Imputation Method. On one hand, the local structure of the empirical data-set105 is preserved since the data are forced to belong to specific manifolds, which depend on the problem nature. On the other hand, the underlying physics of the problem is imposed via supplementary constraints on the data. The imputation procedure is then performed by using an unsupervised learning algorithm that finds the point that minimizes an, in general context-dependent, weighted110 quadratic error, while preserving the local and global physics of the problem. 5 2. Mathematical formulation 2.1. Data-driven simulation-based engineering and sciences (DDSBES): A general framework Our aim is to present a methodology that fits within the context of data-115 driven problems. In particular, problems in which the governing equations of the system are fundamental laws of Physics and the modeling strategy is replaced by data, that, in general, may be incomplete. This means that computations will be carried out directly from data without an a priori parametrization step of the state equations (e.g. constitutive model in Continuum Mechanics).120 In fact, any physical system can be defined as a manifold M, the state space, that corresponds to the admissible states that fulfil a set of equations defining the particular physical problem in hands [39]. Usually, the state space is treated as an embedded manifold M ⊂ Rnin a higher Euclidean dimensional space and is defined in terms of a set of state governing equations F(x) = 0,x∈Rn.125 Observables of the system are magnitudes that are related to the state variables by means of geometric or physical relations. For example, in Continuum Mechanics, forces are observable variables related to the stress components, while displacements are observable variables associated with strains. The relationship between forces and stresses is defined by means of the equilibrium equations,130 while the relation between displacements and strains is derived from kinematic conditions. A set of measurements of the system is a set of observables obtained in particular conditions. For instance, in Continuum Mechanics, the state variables are the stress and the strain tensors with 6 + 6 = 12 components in 3D problems for each135 spatial point (assuming balance of angular momentum yields). But these components are not independent. Indeed, they are related to the observables by kinematic and equilibrium (linear momentum conservation) constraints (note that the standard constitutive relations are implicit in the measurements so they are not considered explicitly).140 The fundamental problem in Physics is stated as: from a set of measure- 6 ments in the physical system given by the values of several observables, derive the value of another observable of interest or/and the rest of state variables. Unfortunately, this inference strongly depends on the quality of measurements as well as on the particular complexity of the system.145 As we have seen, measurements (observables) are related to the state-space variables. However, the information given by these measurements is frequently less informative than the state variables themselves. In other words, measurements are known variables living in a lower dimensional space. For example, the displacements along a given direction do not give us the information about150 the complete strain tensor while the forces over a surface do not characterize the whole stress tensor. In a data-driven framework, the state of the system as well as any desired observable has to be derived from a set of measurements. However, in practical situations, we do not have complete information about the state variables, but only particular measurements in particular states. That155 is, we have an incomplete set of data. As these measurements are related to the state variables, they may be formulated in terms of manifolds. For the Continuum Mechanics problem, for instance, the knowledge of the stress associated with a given plane orientation tells us the relationship between components of the stress tensor, so we have a set of measurements with a reduced dimension of160 the state space. It is essential, therefore, to use these values as points of measure sets being then the goal to properly complete these incomplete measurements to perform data-driven computation. The filling data strategy should take into account therefore the following assumptions: •It should take into account all the measurements, that are formulated165 in terms of manifolds embedded in the state space (incomplete measurements) better than in terms of points (complete measurements). We call this condition the generalization assumption. •It should respect the physics of the system: the new derived data obtained from the incomplete data should be consistent with the geometric struc-170 ture defined by the data manifolds. We call this condition the consistency 7 assumption. •It should guarantee accuracy for the states associated with the input observables used as starting points for predictions. Among all observations in our data-set, those closer to our physical constraints or to our actual ob-175 servables should be overweighted. This is done by means of an appropriate weighting strategy. The methodology for incomplete data processing presented herein should be defined to be in accordance with this framework for DDSBES. First, it is fundamental to define an averaging technique that takes into ac-180 count complete and incomplete measurements according to the generalization assumption and that respects the data structure, according to the consistency assumption. Next, an appropriate weighting strategy for this kind of problems will be defined, i.e, how to compute a set of weighting values in the averaging process in order to make good predictions on the state of the system. Con-185 sequently, the data completion step should take into account how far are the measurements from the known observables and/or the physical manifolds. 8 2.2.3. Filling data using the consistent generalized mean and variance: an unsupervised learning technique The idea of using both generalizations at the same time for missing data295 techniques is natural. Let us assume that we have a set of partially incomplete data Dof size N, that is a set of manifolds {Mj, j = 1,· · · , N}with Mj⊂Rn with respective weights {wj, j = 1,· · · , N}. Here, Rnis the embedding space and Mjrepresents the (incomplete) measurements related to different states. One strategy for data completion of missing data is using the tools presented300 in section 2.2.1. When those measurements are defined in terms of linear manifolds, we have derived, also, a closed linear expression (Proposition 2.1 and Corollary 2.1). The idea is, therefore, to solve the next minimization problem for each measurement Mi,i= 1,· · · , N: min x∈Mi mse(x) = N X j=1 wjd2(x, Mj)(21) The solution xiof this problem is the completed data associated to the305 incomplete data Mi. To summarize, and for linear manifolds, a strategy is assumed for finding a minimum candidate, which the natural and strongly reduces the time required to obtain the actual solution of the minimization problem 21 directly using minimization algorithms. This strategy follows two steps: 1. Global computation step. Computation of the solution point for the310 unconstrained minimization problem using the expression given in Proposition 2.1. 2. Projection step. Computation of the solution point for the constrained minimization problem using the projection of the mean value and uncertainty defined in Corollary 2.1.315 The geometric interpretation of this method is provided in Figure 1. Remark. This two-step method may be used when working with non-linear manifolds as an iterative tangent based algorithm. At each step, the tangent space to the nonlinear manifold is computed at the current point and the linear 15 Figure 1: Geometric representation of the method. A,Band Crepresent three linear manifolds of dimension 2, 1 and 0 respectively, associated to three measurements in a space of dimension 3, the last one complete. Point Pis the generalized mean in the sense defined in Equation 2, that is, the point minimizing the sum of (eventually weighted) squared distances to the manifolds. In orange, the uncertainty ellipsoid, related to the generalized variance, represents the spread of (eventually incomplete) measurements. Points PA,PBand PCare the consistent generalized means associated to each of the manifolds in the sense defined in Equation 21, that is, the projection of the generalized mean on each manifold. In green, the associated uncertainty ellipsoid, related to the consistent generalized variance, for each manifold is depicted. 16 problem is solved. Thus, a point belonging to the tangent space is obtained.320 Using the exponential map [41], it is possible to obtain an associated point belonging to the manifold, which is used as the starting point for the next iteration. This strategy is usual in nonlinear computational mechanics [42, 43]. The iteration scheme stops when the distance between the subsequent global solutions is lower than a given tolerance. This construction is illustrated by325 the schematic diagram in Figure 2. This algorithm may be computationally expensive and is very dependent on the manifold smoothness and convexity, being this type of problems out of the scope of this paper.  In any case, when the algorithm achieves convergence, Nvalues of Xare obtained, one for each manifold, xi,i= 1,· · · , N, representing the expected330 value associated with the manifold Mi. Besides, for linear problems, we obtain for each manifold a value s2 Mithat characterizes the uncertainty related to the manifold Miand the matrix SMicharacterizing the uncertainty spread on this manifold. This uncertainty ellipsoid can be seen in Figure 1 for linear manifolds of dimension 0, 1 and 2, respectively.335 17 Figure 2: Extension of the method to nonlinear problems: Initially, for each manifold (A,B and C) an initial point is selected (A0,B0and C0=C) and the tangent spaces are computed. Then, the optimal completion, i.e., the generalized consistent mean for each linear manifold (in red), is computed. Using the exponential map, these projections on the linear manifolds are translated to the respective manifolds (A,Band C) obtaining a new point in each manifold (A1,B1and C1=C0). The process is repeated until convergence. 18 2.3. Introducing physical laws: weighting strategy Here we introduce the weighting strategy. Let us consider a system defined in terms of a physical manifold Mand a set of measure manifolds Nj, j = 1,· · · , N. The starting point from which we want to derive the state of the system is another measure, that is, another manifold Nthat could be, for340 instance, boundary conditions for a given problem. The fundamental idea is to compute the manifolds Mj=Nj∩ N ∩ M =Nj∩ P where P=N ∩ M and to perform the unsupervised learning strategy for this reduced space. This strategy has two direct consequences: •Physical consequence: Since we are only learning second order statistics345 of physically admissible manifolds, the result will have a more physical sense. •Numerical consequence: Projections are performed in a smaller space so the computational cost will be lower. In that case, we are looking for a physically admissible (incomplete) data350 point measure that has the lowest uncertainty. However, this strategy can dramatically fail for few data with non-negligible uncertainty. For example, a measure manifold Njmay be close to the manifold Pbut P ∩Nj=∅so this measure will not be used for the system learning, even though it is very close to the real state.355 An intermediate solution is using an activation function in the learning step, depending on the distance to the manifold P,d=d(P,Nj). That is, a function φ:R+→[0; 1] so that if z= 0, φ(z) = 1 and eventually if z→+∞, φ(z)→0. Given Nj⊂Rn, defining uj=u(Nj) = d(Nj,P), it is possible to define wj=φ(uj) in the learning process. In other words, the nearer the con-360 sidered data set to P, the higher the weight should be in the minimization of the optimal distance. Some possible activation functions are the step function φ(u) = χ[0;a](u), where χAis the characteristic function of the set Aand a≥0, radial basis functions (RBF), φ(u) = exp −u2 2ς2with ς > 0, homographic 19 functions, φ(u) = a a+bu k, with a, b, k > 0, or generalized ramp functions365 φ(u) = 1−u akχ[0;a](u), with a, k > 0. Figure 3 illustrates the geometric idea under the presented filling method when combined with the physics of the problem. 20 Figure 3: Geometric idea of the learning phase for data completion of the incomplete measure N1and its corresponding uncertainty. Blue manifolds Nj,j= 1,··· ,5.represent measurements (all incomplete except N5, which is a point). These measurements may correspond to states far from the current state of the system, that should belong to the orange physical manifold P=M ∩ N that corresponds to the intersection of the points satisfying the governing equations of the problem (equilibrium, thermodynamics, Maxwell’s equations etc.) and knowledge about the system state, e.g. boundary conditions or measured control variables. Each of the measure manifolds is at a certain distance ujfrom the physical manifold. In particular, in the figure case, u1= 0 and u3= 0 because N1∩ M 6=∅and N3∩ M 6=∅. Note that this situation is frequent in a three-dimensional state-space, but it is less and less probable when the dimension of the total space becomes very high. From these distances, the weights are calculated by wj=φ(uj), so that if u= 0, φ(u) = 1 and if u→ ∞,φ(u)→0. Thus, if x∈R3, the distance D2(x) = P5 j=1 wjd2 j(x) = P5 j=1 φ(uj)d2(x, Nj) is minimized (generalized mean) and the solution point P(in red) is projected onto the measure manifolds (consistent mean) obtaining the measure completion, as illustrated in Figure 1. Likewise, depending on the ”spread” of the sets, an ellipsoid of three-dimensional uncertainty (in red), is obtained, related to the generalized variance. This ellipsoid is projected in each of the manifolds, obtaining ellipsoids related to the associated consistent variance (in green). Note that if the red ellipsoid is very slender in the direction orthogonal to the measure manifold that we are completing, this would have no impact on the projected green ellipsoid. The point P1, is then the generalized consistent mean associated to incomplete measure N1weighted by the neighborhood to the problem physics. 21 2.4. Looking for the nearest measure to a given point: solving the data-driven problem370 2.4.1. Preliminary mathematical results Let M ⊂ Rnbe an embedded manifold of the Euclidean space of dimension nwith associated probability distribution ρ, that accounts for the probability distribution of a given random point xbelonging to M. For instance, in the case of linear manifolds, Mcan be described as M=p+hv1,v2,· · · ,vmi,375 where hv1,v2,· · · ,vmiis the linear span generated by v1,· · · ,vm. That is R={p;v1,v2,· · · ,vm}is a basis of Mand m= dim(M). Let x∈ M and ρ:M → R+, be the probability distribution describing the position of x. Then we define the square distance random variable D2= d2(p, X), with p(deterministic) and X(random) defined by their coordinates380 p,Xin a (global) reference frame, as D2=||p−X||2(22) It is possible to define (under some integrability conditions) the expected value E[D2] = RMd2(p, X)ρ dV , or, using the point coordinates x(note that Xis a random vector while xthe vector point coordinates), E[D2] = RM||p− x||2ρ(x)dV (x). More important than the explicit computation of the expected385 value of D2in terms of a given parametrization are the following results for linear manifolds, proven in Appendix A, and relating the moments of D2to the moments of X. As usual, the random variable Xand the point pare identified with their coordinates description X, and p. If µis the expected value of the random vector X(µi=E[Xi]), Σits variance-covariance matrix390 (Σij =E[XiXj]−E[Xi]E[Xj]) and Υits fourth order moment tensor (Υijkl = E[XiXjXkXl]), then: E[D2] = ||p−πM(p)||2+||πM(p)−µ||2+ Tr(Σ) (23) Var[D2] = I:Υ:I−Tr(Σ)+(πM(p)−µ)T(πM(p)−µ)(24) 22 Here, Iis the second order identity tensor. Finally, under normality conditions, D2follows a noncentral χ2distribution with m= dim(M) degrees of freedom and non-centrality parameter λ= (πM(p)−µ)T(Σ)−1(πM(p)−µ).395 2.4.2. The closest point in a stochastic sense Let us go back now to the methodology and tools introduced in section 2.2. Once the constrained minimization problem is solved for each incomplete measurement Miand all filled data points xiare derived, it is possible to compute how far is a state p∈Rnof the system from a given data-point i.400 Moreover, we can define for each p∈Rnwhich is the closest measure, and from this to define a tessellation of the state space in terms of the measurements. One could consider the deterministic distance di=d(p, xi) but this distance would not have into consideration the accuracy of the filling step and the effect on physical weights on uncertainty. It is more natural to consider a stochastic405 distance. Indeed, considering again the random variable D2 idefined at Section 2.4.1 associated with the manifold Mi, we may define, denoting s2 Miby s2 iand πMiby πi: d2 i=E[D2 i] = d2(p, πi(p)) + d2(πi(p), xi) + s2 i(25) The manifold Miverifying that d2 iis minimal is the closest manifold in the statistical sense to the point p. Besides, each term in d2 ihas its own410 interpretation: •Ti,1=d2(p, πi(p)) is the statistical error due to finite measurements of the sample. It is related to the lack of knowledge about the system, since the information is obtained by means of a given finite data-set. The more measurements are added, the lower the error usually is.415 •Ti,2=d2(πi(p), xi) is inherent to the manifold and depends on the manifold selection. It is unavoidable to some extent. •Ti,3=s2 iis the term associated with the uncertainty and is characteristic of the self-learning process: the worse the manifold learning, the higher 23 this term. Locating properly the real point in a manifold, even though it420 actually belongs to that manifold, is less accurate when this term increases. Moreover, if measurement uncertainty is taken into consideration, it is possible to state Ti,3=s2 i+s02 iwhere s02 iis the quadratic uncertainty of the i-th measurement and, therefore, orthogonal to the uncertainty associated to the filling procedure. This uncertainty, nevertheless, is not being considered in the425 applications presented in this work. The geometric idea behind these considerations is illustrated in Figure 4. Once the learning step is finished, the data-driven problem may be solved as usual [15, 17], provided the uncertainty of the completion is considered, as just explained. Figure 5 illustrates the geometric idea.430 24 methods. Note that the error of the presented method (UL) does not increase with the amount of lost data. This is due to the fact that the performance of the method depends on how far are the missed data points and not on their number, i.e. data quality and not data quantity. On the other hand, standard interpolation techniques are rather dependent on the distance of the missing510 data from the true solution than in the volume of missing data itself. Moreover, low order interpolation techniques (nearest or linear interpolation) are sometimes unable to reproduce the underlying data structure, whereas high order interpolation techniques (cubic splines) performance is strongly dependent on data sampling [44]. Shape preserving interpolation, for example,515 was conceived as a compromise solution to these problems, but is still strongly dependent on the missing data process as has been demonstrated [45]. Listwise deletion is the most robust method with respect to missing data fraction, but has statistical power and the bias problems for MNAR data as reported in literature [27], [23]. The proposed method shows a more robust behavior with520 respect to the missing data fraction. 31 3.2. Model-based data-driven problem The performance of the method is now illustrated in a classical problem of strength of materials. Let us consider a two-end clamped beam of length Lunder bending by a linearly distributed load q=q(x), 0 ≤x≤L. In the525 Euler-Bernoulli framework, and supposing the beam composed of a linear elastic material with Young modulus E=E(x), 0 ≤x≤L, and a section with moment of inertia I=I(x), the vertical beam displacement u=u(x), 0 ≤x≤L, may be computed solving the linear differential equation: d2 dx2EI d2u dx2=q(30) with boundary conditions:530 u(0) = du dxx=0 = 0 u(L) = du dxx=L = 0 (31) Equation (30) with boundary conditions (31) may be solved numerically using any standard numerical procedure (e.g. Finite Elements or Finite Differences). Once the problem is discretized using a mesh of characteristic size h and the boundary conditions are applied, the nodal displacements uhare obtained solving the linear system (it is important to note that this equation is535 characteristic of a broad family of linear discretized problems in Physics and Engineering, not only the Euler-Bernoulli bending beam): Khuh=fh(32) In order to test our method, we may proceed by considering the equation (32) as a physical constraint to a data-set of measurements E={uh0 i, i = 1,· · · , N} equally spaced h0. Note that this approach makes sense when h0<< h (this540 may be the case when equation (32) is computationally expensive to solve with 32 very fine meshes while measurements are easy to obtain). In this case, we deal with complete measurements, but they are subjected to error (bias and noise for example due to experimental reasons). The presented method is able to detect how far a given measure is from the physics of the problem in terms of the545 distance to the manifold defined by equation (32). The different measurements will be weighted differently depending on their distance to the manifold. Recall that a standard procedure of averaging all measurements may induce an error if there is a systematic bias in the measurements, which is a well-known problem of mean imputation [25]. The weighting strategy considered associates the bias550 in the estimation to data quality in a physical sense. In order to illustrate the application of the methodology, let us solve the defined problem for q(x) = 10 kN/m, L= 10 m, E(x)I(x)=1·106kN ·m2. In that case, the analytical solution is given by u(x) = −qL4 24EI x L−x2 L22. The measurements are randomly generated from the analytical solution sam-555 pled in a mesh of size h0=1 m−1,uk∼ N(u(x=kL h0) + b, σ), where b=α¯uis a bias, α∼ N(0.03,0.03), σ=β¯u,β= 0.002 and ¯u=1 30 qL4 EI is the mean value of the analytical solution. For this example, m= 100 points are considered along the beam length and N= 6 samples are evaluated. The six families of points are shown in Figure 7. Also, and in the same figure, the results computed560 for different spacings hused to establish the physical constraint are shown. In particular, the mean of the samples is represented by the bold orange color line and the true analytical solution is represented by the continuous bold blue line. The RBF function φ(u) = exp −u2 2ς2with ς= 5 ·10−4m was selected as the weighting function.565 As pointed out above, the method here described corrects partially the bias of the measurements. This correction is done automatically by computing the distance of each measure sample to the physical manifold defined by the discretized equilibrium equation and boundary conditions, and transforming this distance using the RBF function. A finer mesh in the discretized physical prob-570 lem (lower h) takes into account more points to compute the distance from the measure to the physics of the problem, while this physics is more accurate. It 33 012345678910 -3 -2 -1 0 1 2 10-4 Figure 7: Comparison of the solutions obtained using the presented method for different mesh sizes hto define the physical constraint. 34 is clear, however, that the accuracy of the estimation in cases as this one with a systematic positive bias will never be better than the best measure. In order to have a deeper understanding of the method, it is worth to make575 a physical interpretation of this example. For a fixed h, equation (32) is a constraint relating n= 1/h + 1 variables of the m-dimensional space. The manifold Ph≃Rnof dimension nmdefined by these ncoordinates is the manifold where the relevant physics of the problem is evaluated. It is easy to figure up what the method does by considering only the projections of the points580 in Rmon Ph. In Ph, the physical manifold (that is, the constraint) is given by a single point p: the numerical solution obtained solving the equation (32). Therefore, the distances of the different samples to the physical manifold may be interpreted as Euclidean distances in Ph≃Rnbetween sample projections on Ph(i.e. the consideration of the ncoordinates related to the mesh with585 hspacing) and point P,di=||πPh(ui)−uh||.We could have defined the distance in a more general framework, such as Hilbert spaces, using the Finite Element approximation, but a simpler and more interpretable norm was selected for illustration purposes. The weighted mean is then computed by using the solution uhas the reference point in the weights computations, wi=φ(di).590 The physics of the problem is inferred outside the manifold Ph, that is, in (Ph)⊥. In Figure 8 all described geometric elements are shown for h= 1/10: projections over manifold Phof m-dimensional samples, numerical solution uh and associated generalized mean and variance. The analytical solution is also plotted even though it does not appear in the computation.595 35 012345678910 -4 -3 -2 -1 0 1 210-4 Numerical solution Sample 1 Sample 2 Sample 3 Sample 4 Sample 5 Sample 6 Analytical sol. Gen. mean and std. dev. band Figure 8: Projections over the manifold Phof all geometric elements used in computations for h= 1/10. The dotted magenta line represents the generalized mean (in this case it is equivalent to the weighted mean with complete measurements). The error bands correspond to the projections over the different coordinates lines (manifolds of dimension 1) of the generalized variance ellipsoid (in orange in Figure 1 or Figure 3). Green ellipsoids represented in these figures collapse in the data samples points because samples are complete so the measure manifolds are points. 36 3.3. Model-free data-driven problem The final application lays within the framework in which the method has been conceived: model-free data-driven problems. The objective is now to solve a real physically-based problem, formulated in terms of a set of governing equations encoding the physical information, and some empirical knowledge, formu-600 lated in terms of a given data-set. The main objective of this example, despite its simplicity, is to illustrate the performance and some of the properties and capabilities of the technique proposed in DDSBES rather than comparing it with other filling data methods, as in the previous examples. With this aim, the problem is first formulated in605 a classical framework approach, where no data-set is considered and the entire physics of the problem is supposed to be known and parametrized. Secondly, the problem is reformulated in the data-driven framework, where only the sound physics is postulated and the rest of the physical structure is built from data. Let us suppose two reservoirs connected by a channel with section S, hy-610 draulic diameter Dand length Las illustrated in Figure 9. A fluid flows from one reservoir to the other depending on the water level in each of them, y1for the reservoir 1 and y2for the reservoir 2. The section of each reservoir is defined for each height yiby a function Si=Si(yi), i= 1,2, being clear that the volume occupied by the fluid at reservoir ifor a height yiis Vi(yi) = Ryi 0Si(u)du.615 When considering the physics of the problem, two sound laws are invoked: conservation of mass (1) and conservation of energy (2). The first one is equivalent to impose (under the assumption of fluid incompressibility) zero net flow, Qnet = 0. That writes: dV1 dt +dV2 dt = 0 (33) Using the expression of Viin terms of yiwe get:620 S1(y1) ˙y1+S2(y2) ˙y2= 0 (34) The second considers the energy loss due to viscous dissipation, unless we 37 Figure 9: Schematic diagram of the fluidic device used to validate the presented methodology. consider an inviscid fluid. Using Bernoulli energy conservation statement in terms of water height, conservation of energy writes y1−y2=IL (35) where Iis the hydraulic head loss slope, I=I(Q) = fD1 2g Q2 S2Dand Q= S1(y1) ˙y1. In order to close the system of equations (34) and (35), it is necessary625 to define the Darcy friction factor fDor at least to express it in terms of the state space variables y1,y2, ˙y1and ˙y2, which is the critical step in classical approaches. In general, either one additional hypothesis is assumed (for example laminar regime), carry out simulations with complex fluid flow models or use a semi-empirical equation such as K´arm´an -Prandtl resistance equation for the630 smooth turbulent regime, [46], [47], the well-known Colebrook-White equation [48] or other more recent equations for the transition from a smooth pipe to a rough pipe flow [49], [50], [51]. One way or another, these approaches complete the physics based on the particular hypothesis stated a priori. Assuming a laminar regime, fD=64 Re =64Sν QD , with νthe kinematic viscosity635 and defining the initials conditions, the problem may be solved numerically integrating the system of equations. For instance, we may fix y1(t= 0) = H, y2(t= 0) = 0, ˙y1(t= 0) = 0 and ˙y2(t= 0) = 0 (reservoir 1 at level Hand 38 reservoir 2 empty, both at rest). The data-driven approach is based on the use of a data-set that will com-640 plete implicitly the physics of the problem a posteriori. Using this approach, only the conservation of mass (equation (34)) is taken into account while the energy equation is replaced by a data-set sampled from the state-space M= {(y1, y2,˙y1,˙y2)|S1(y1) ˙y1+S2(y2) ˙y2= 0}. We have therefore a data-set S= {(yi 1, yi 2,˙y1i,˙y2i)|i= 1,· · · , N}. As we deal again with a missing data problem,645 we define the missing process as follows: due to experimental limitations, it is impossible to measure the velocity of the free surface level when the free surface is higher than a defined threshold h∗. This condition tries to reproduce a realistic missing data process related to the experimental setup. We are clearly in a MNAR situation. The presented method, which is local in a certain sense650 due to the weighting process, should be insensitive to the missing process. In order to illustrate the methodology, let us solve the problem in a particular case using both the classical approach assuming laminar regime and the datadriven framework, using as data-set the solution in the laminar regime with noise: for each variable xnoise =xexact +where ∼ N (0, pσ), being p= 1/10655 and σthe standard deviation of the exact values xexact. A truncated cone geometry is assumed for both reservoirs with lower radii r= 1.33 mm and slope s= 0.035. The density and dynamic viscosity of water are ρ= 1000 kg/m3 and µ= 1.006 ·10−3Pa ·s, ν=µ ρand g= 9.81 m/s2, respectively. The initial level of the fluid is H= 5 cm and the channel has length L= 10 cm and a660 rectangular section of width w= 750 µm and height h= 200 µm. The laminar solution in terms of the two heights and the flow is shown in Figure 10. Our filling data strategy is tested by solving the data-driven problem with a loss fraction of data defined by fixing the threshold h∗as described before. The components y1and y2of the complete data-set are shown in Figure 11a.665 Our aim is solving the problem when we observe a value of y1=m, that is, to obtain the other state variables y2, ˙y1and ˙y2. In other words, to find the closest point to our incomplete data-set satisfying mass conservation (physical manifold M) and using as reference our initial set of observations (observation 39 0 100 200 300 400 500 600 700 800 900 1000 0 10 20 30 40 50 Channel of width 0.75 mm and height 0.20 mm y1 y2 0 100 200 300 400 500 600 700 800 900 1000 0 50 100 150 Figure 10: Classical solution of the problem solving the system of ordinary differential equations. manifold N). The manifold Mis defined by the mass conservation equation670 (34) while the manifold Nis defined by equation y1−m= 0. Incomplete data are filled using the methodology presented in this work. Completed data are shown in Figure 11b for m= 42 mm and h∗= 1.1·m. As stated before, an RBF function is chosen as weighting function with ς= 0.2 mm while, for illustration purposes, weights are computed only in terms of the distance to the675 manifold Ninstead of N ∩ M, which is nonlinear. Finally, the solution point is chosen by solving the data-driven minimization problem suggested in [16, 17], that is, by minimizing the distance to the (filled) data-set, provided the physical constraint. The difference here is the fact that a stochastic distance is computed, including the deterministic term (Ti,1and Ti,2, related to the generalized mean)680 and the quadratic uncertainty term (Ti,3, related to the generalized variance), both resulting from the filling procedure, as it is described in section 2.4.2. Figure 11b depicts the solution of the problem in terms of intuitive plots or state variables but does not illustrate the geometry behind. In order to 40 Appendix A. Mathematical proofs Proposition Appendix A.1 (Expected value of the squared distance). Let Va linear manifold and let πVthe orthogonal projection on V,P∈Rn, X∈ V a random vector and D=d(P, X). We identify Xwith random vector795 Xand point Pwith its coordinates pin a given reference frame. Let us suppose that µ=E[X]and Σ= COV(X)are finite. Therefore: E[D2] = ||p−πV(p)||2+||πV(p)−µ||2+ Tr(Σ) (Appendix A.1) Proof:. Using Pythagoras theorem we have: D2=||p−X||2 =||p−πV(p)||2+||πV(p)−X||2(Appendix A.2) Therefore E[D2] = E||p−πV(p)||2+||πV(p)−X||2 =||p−πV(p)||2+E||πV(p)−X||2(Appendix A.3) Now, we have800 E||πV(p)−X||2= Tr(Σ) + ||E[πV(p)−X]||2 = Tr(Σ) + ||πV(p)−E[X]||2(Appendix A.4) Combining equations Appendix A.3 and Appendix A.4 we obtain the result.  We may observe that when using a reference frame in V, rang(Σ) = dim(V) and then Tr(Σ) has as many terms as the dimension of V. An analogous result for the variance of D2can be derived.805 47 Proposition Appendix A.2 (Variance of the squared distance). Let Va linear manifold and let πVthe orthogonal projection on V,P∈Rn,X∈ V a random vector and D=d(P, X). We identify Xwith random vector Xand point Pwith its coordinates pin a given reference frame. Let us suppose that µ=E[X],Σ= COV(X)and Υ=Υ(X)the centered tensor moment of order810 4 of Xare finite. Therefore: Var[D2] = I:Υ:I−hTr(Σ)+(πV(p)−µ)T(πV(p)−µ)i(Appendix A.5) Proof:. We have seen in previous proof that: D2=||p−πV(p)||2+||πV(p)−X||2(Appendix A.6) Therefore: Var[D2] = Var ||p−πV(p)||2+||πV(p)−X||2 = Var ||πV(p)−X||2(Appendix A.7) But we have, Var ||πV(p)−X||2=I:Υ:I−Tr(Σ) + E[πV(p)−X]TE[πV(p)−X] =I:Υ:I−Tr(Σ)+(πV(p)−µ)T(πV(p)−µ) (Appendix A.8) Combining equations Appendix A.7 and Appendix A.8 we obtain the result.815  Under normality conditions, we have the following result: Proposition Appendix A.3 (Squared distance distributional properties). Let Va linear manifold and let πVthe orthogonal projection on V,P∈Rn, X∈ V a random vector and D=d(P, X). We identify Xwith random vector820 48 Xand point Pwith its coordinates pin a given reference frame. Let us assume that Xfollows a multivariate normal distribution with mean µ=E[X]and covariance matrix Σ= COV(X). Let χ2= (p−X)T(Σ)−1(p−X). Then χ2 follows a non-central χ2distribution with k= dim(V)degrees of freedom and non-centrality parameter λ= (πV(p)−µ)T(Σ)−1(πV(p)−µ).825 Proof:. Let U=p−X. Therefore, Ufollow a multivariate normal distribution with mean p−µand covariance matrix Σ. Then, Σ−1/2U∼ N(Σ−1/2(p−µ),I). Using non-central χ2distribution definition: χ2=UTΣ−1U= (Σ−1/2U)T(Σ−1/2U)∼χ2(n, λ) (Appendix A.9) Where n= rang(Σ) = dim(V) and λ= (Σ−1/2(p−µ))T(Σ−1/2(p−µ)) = (p−µ)TΣ−1(p−µ)).830 But we have, p−µ= (p−πV(p)) + (πV(p)) −µ). As p−πV(p)∈ ker(Σ−1),then (p−µ)TΣ−1(p−µ)=(πV(p)−µ)TΣ−1(πV(p)−µ).  49 References [1] L. Atzori, A. Iera, G. Morabito, The internet of things: A survey, Computer835 networks 54 (15) (2010) 2787–2805. [2] J. Manyika, M. Chui, B. Brown, J. Bughin, R. Dobbs, C. Roxburgh, A. H. Byers, Big data: The next frontier for innovation, competition, and productivity. [3] D. T. Larose, Discovering knowledge in data: an introduction to data min-840 ing, John Wiley & Sons, 2014. [4] T. M. Mitchell, Machine learning. 1997, Burr Ridge, IL: McGraw Hill 45 (37) (1997) 870–877. [5] S. Khan, Introduction to machine learning (adaptive computation and machine learning series), Natural Language Engineering 14 (01) (2008) 133–845 137. [6] V. Vapnik, The nature of statistical learning theory, Springer science & business media, 2013. [7] F. Rosenblatt, The perceptron: A probabilistic model for information storage and organization in the brain., Psychological review 65 (6) (1958) 386.850 [8] A. Krenker, J. Bester, A. Kos, Introduction to the artificial neural networks, in: Artificial neural networks-methodological advances and biomedical applications, InTech, 2011. [9] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, et al.,855 Mastering the game of go with deep neural networks and tree search, Nature 529 (7587) (2016) 484–489. [10] F. Darema, Dynamic data driven applications systems: A new paradigm for application simulations and measurements, in: International Conference on Computational Science, Springer, 2004, pp. 662–669.860 50 [11] R. E. Kalman, A new approach to linear filtering and prediction problems, Journal of basic Engineering 82 (1) (1960) 35–45. [12] B. Peherstorfer, K. Willcox, Dynamic data-driven reduced-order models, Computer Methods in Applied Mechanics and Engineering 291 (2015) 21– 41.865 [13] R. Ibanez, E. Abisset-Chavanne, J. V. Aguado, D. Gonzalez, E. Cueto, F. Chinesta, A manifold learning approach to data-driven computational elasticity and inelasticity, Archives of Computational Methods in Engineering (2016) 1–11. [14] P. Ladev`eze, The large time increment method for the analysis of struc-870 tures with non-linear behavior described by internal variables, COMPTES RENDUS DE L ACADEMIE DES SCIENCES SERIE II 309 (11) (1989) 1095–1099. [15] T. Kirchdoerfer, M. Ortiz, Data-driven computational mechanics, Computer Methods in Applied Mechanics and Engineering 304 (2016) 81–101.875 [16] T. Kirchdoerfer, M. Ortiz, Data driven computing with noisy material data sets, arXiv preprint arXiv:1702.01574. [17] 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, Computer Methods in Applied Mechanics and Engi-880 neering 328 (2018) 752–774. [18] D. B. Rubin, Inference and missing data, Biometrika 63 (3) (1976) 581–592. [19] J. L. Schafer, J. W. Graham, Missing data: our view of the state of the art., Psychological methods 7 (2) (2002) 147. [20] P. D. Allison, Missing data: Quantitative applications in the social sciences,885 British Journal of Mathematical and Statistical Psychology 55 (1) (2002) 193–196. 51 [21] P. D. Allison, Missing data techniques for structural equation modeling., Journal of abnormal psychology 112 (4) (2003) 545. [22] J. W. Graham, Missing data analysis: Making it work in the real world,890 Annual review of psychology 60 (2009) 549–576. [23] M. Nakai, W. Ke, Review of the methods for handling missing data in longitudinal data analysis, International Journal of Mathematical Analysis 5 (1) (2011) 1–13. [24] C.-Y. J. Peng, M. Harwell, S.-M. Liou, L. H. Ehman, et al., Advances in895 missing data methods and implications for educational research, Real data analysis 3178. [25] J. L. Schafer, Analysis of incomplete multivariate data, CRC press, 1997. [26] P. D. Allison, Missing data: Sage university papers series on quantitative applications in the social sciences (07–136), Thousand Oaks, CA.900 [27] A. Briggs, T. Clark, J. Wolstenholme, P. Clarke, Missing.... presumed at random: cost-analysis of incomplete data, Health economics 12 (5) (2003) 377–392. [28] H. Junninen, H. Niska, K. Tuppurainen, J. Ruuskanen, M. Kolehmainen, Methods for imputation of missing values in air quality data sets, Atmo-905 spheric Environment 38 (18) (2004) 2895–2907. [29] R. J. Little, D. B. Rubin, Bayes and multiple imputation, Statistical Analysis with Missing Data, Second Edition (2002) 200–220. [30] D. B. Rubin, Multiple imputation after 18+ years, Journal of the American statistical Association 91 (434) (1996) 473–489.910 [31] N. Schenker, J. M. Taylor, Partially parametric techniques for multiple imputation, Computational statistics & data analysis 22 (4) (1996) 425– 446. 52 [32] J. L. Schafer, Multiple imputation: a primer, Statistical methods in medical research 8 (1) (1999) 3–15.915 [33] H. Hartley, R. Hocking, The analysis of incomplete data, Biometrics (1971) 783–823. [34] C. K. Enders, A primer on maximum likelihood algorithms available for use with missing data, Structural Equation Modeling 8 (1) (2001) 128–141. [35] C. K. Enders, D. L. Bandalos, The relative performance of full informa-920 tion maximum likelihood estimation for missing data in structural equation models, Structural equation modeling 8 (3) (2001) 430–457. [36] A. P. Dempster, N. M. Laird, D. B. Rubin, Maximum likelihood from incomplete data via the em algorithm, Journal of the royal statistical society. Series B (methodological) (1977) 1–38.925 [37] C. K. Enders, Applied missing data analysis, Guilford Press, 2010. [38] B. G. Tabachnick, L. S. Fidell, S. J. Osterlind, Using multivariate statistics. [39] G. J. Sussman, J. Wisdom, Structure and interpretation of classical mechanics, Mit Press, 2015. [40] A. Papoulis, Probability & statistics, Vol. 2, Prentice-Hall Englewood930 Cliffs, 1990. [41] M. P. Do Carmo, J. Flaherty Francis, Riemannian geometry, Vol. 115, Birkh¨auser Boston, 1992. [42] J. Simo, N. Tarnow, M. Doblare, Non-linear dynamics of three-dimensional rods: Exact energy and momentum conserving algorithms, International935 Journal for Numerical Methods in Engineering 38 (9) (1995) 1431–1473. [43] S. Niroomandi, I. Alfaro Ruiz, E. Cueto Prendes, Real-time simulation of surgery by model reduction and x-fem techniques. 53 [44] D. Kahaner, C. Moler, S. Nash, Numerical methods and software, Englewood Cliffs: Prentice Hall, 1989.940 [45] F. N. Fritsch, R. E. Carlson, Monotone piecewise cubic interpolation, SIAM Journal on Numerical Analysis 17 (2) (1980) 238–246. [46] H. J. Tracy, C. Lester, Resistance coefficients and velocity distribution smooth rectangular channel, US Government Printing Office, 1961. [47] B. J. McKEON, M. V. ZAGAROLA, A. J. SMITS, A new friction factor945 relationship for fully developed pipe flow, Journal of Fluid Mechanics 538 (2005) 429443. doi:10.1017/S0022112005005501. [48] C. F. Colebrook, T. Blench, H. Chatley, E. Essex, J. Finniecome, G. Lacey, J. Williamson, G. Macdonald, Correspondence. turbulent flow in pipes, with particular reference to the transition region between the smooth and950 rough pipe laws.(includes plates)., Journal of the Institution of Civil engineers 12 (8) (1939) 393–422. [49] N. Afzal, Friction factor directly from transitional roughness in a turbulent pipe flow, Journal of Fluids Engineering 129 (10) (2007) 1255–1267. [50] N. Afzal, Erratum:” friction factor directly from transitional roughness in955 a turbulent pipe flow”[asme trans. j. fluid eng., 2007, 129, pp. 1255-1267], Transactions of the ASME-I-Journal of Fluids Engineering 133 (10) (2011) 107001. [51] N. Afzal, A. Seena, A. Bushra, Turbulent flow in a machine honed rough pipe for large reynolds numbers: General roughness scaling laws, Journal960 of hydro-environment research 7 (1) (2013) 81–90. 54