Simulation of any physical process requires definition of the physical model, method – analytical or numerical, to solve the set of equations describing the physical model and the parameters expressing the body properties and boundary conditions. This paper focus on two latter aspects of the numerical simulation process. Precise determination of the model quantities are crucial for high quality of the model predictions and accurate reflection of real system. Determination of the process parameters is defined as an inverse problem. Following this the sensitivity analysis is applied as the preliminary step of the inverse analysis to reduce the number of model evaluations and to increase the inverse calculations robustness and efficiency. Sensitivity analysis techniques show how "sensitive" is a model to its input parameters variations and to changes of the model structure. As the example the sensitivity analysis was applied to the 2D DC borehole resistivity measurements simulation problem solved with hp-Finite Element Method

Szeliga, D.

English

UPCommons. Portal del coneixement obert de la UPC

1359   APPLICATION OF SENSITIVITY ANALYSIS – PRELIMINARY STEP OF THE PROCESS PARAMETERS ESTIMATION D. SZELIGA Department of Applied Computer Science and Modelling Faculty of Metals Engineering and Industrial Computer Science Akademia Górniczo-Hutnicza Al. Mickiewicza 30 30-059 Kraków, Poland e-mail: szeliga@agh.edu.pl Key words: Sensitivity analysis, inverse problem, parameter identification, hp-adaptivity, electrical resistivity Abstract. Simulation of any physical process requires definition of the physical model, method – analytical or numerical, to solve the set of equations describing the physical model and the parameters expressing the body properties and boundary conditions. This paper focus on two latter aspects of the numerical simulation process. Precise determination of the model quantities are crucial for high quality of the model predictions and accurate reflection of real system. Determination of the process parameters is defined as an inverse problem. Following this the sensitivity analysis is applied as the preliminary step of the inverse analysis to reduce the number of model evaluations and to increase the inverse calculations robustness and efficiency. Sensitivity analysis techniques show how "sensitive" is a model to its input parameters variations and to changes of the model structure. As the example the sensitivity analysis was applied to the 2D DC borehole resistivity measurements simulation problem solved with hp-Finite Element Method. 1 INTRODUCTION Modeling of any physical problem requires precise quantitative information of the model parameters. Some of them are derived from physical laws, others are of phenomenological nature. Proper physical and mathematical description of the problem as well as selection of the solution method and accurate estimation of the model parameters are crucial for the high quality of the modeling results. The paper focuses on the problem of model parameters estimation and the efficient methods to determine the parameters.  Most problems describing physical phenomena related to identification of some quantities are defined as inverse problems [1],[2]. Those problems are hard to solve due to non-unique solution and the lack of the model output stability with respect the identified parameters. Another aspect is efficiency of the identification. Models of physical phenomena are based on differential equations and solved with time consuming numerical methods (e.g. finite element method, finite volume method, particles method). All those features motivate to develop the robust parameter identification method of high efficiency with respect to the calculation time. Classical inverse method was developed and applied by the Author to identify rheological XI International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XI E. Oñate, D.R.J. Owen, D. Peric and B. Suárez (Eds) 1360D. SZELIGA  2material properties [3]–[5]. One of the main disadvantages of that approach is calculation time. The idea of the modified inverse method is to supply the classical algorithm with sensitivity analysis as the preliminary step of the solution to decrease the calculation time. As the application example problem of DC (direct current) borehole resistivity measurements is analyzed. 2 PROBLEM FORMULATION 2.1 Direct and inverse problem Integral or differential equations describing any physical phenomena are set out in terms of functional analysis as: :K X Y→  (1) where X and Y are normed spaces and K is a mapping (linear or nonlinear).  The direct problem is formulated as evaluating y = K(x) ∈ Y for given x ∈ X and an operator K that is equivalent to solve a boundary value problem for differential equation or to evaluate an integral. The inverse problem is defined as evaluating the x ∈ X value for given K and y ∈ Y [1],[2]. It could be shown that inverse problems described as the integral/differential equations are ill-posed in the sense of Hadamard [6]. Those problems require regularization procedure and one of the solution is transforming them to the following, well-posed, problems: 2x Kx y−  (2) The form (2) leads to minimization with respect to the parameters, which are identified: boundary conditions parameters or material parameters. In terms of optimization terminology, the inverse problem is to find the minimum of the objective function: ( ) 2:x Kx yδΦ = −  (3) where yδ is the perturbated (measured) data such that || ||y yδ δ− ≤ , y ∈ K(X) – exact solution of equation (1).  The objective function (3) depends on the norm in the Y space, if it is supplied with Euclid’s norm, the objective function is defined as an average square root error (the Euclidean distance) between calculated and measured quantities. Equation (1) can be expressed as: ( )* *x I aK K x aK y= − +  (4) where K* is adjoint operator and a>0 is the a number. For equation (4) the iteration procedure scheme : ( )0* 1 *: 0, 1, 2,...m mxx I aK K x aK y m−== − + = (5) 1361D. SZELIGA  3leads to the solution of the inverse problem described by (2). The general flow chart of the inverse analysis algorithm is presented in figure 1. Regardless of the process type inverse problem and independently of the applied method of the solution, the algorithm consists of three parts: • The set of process outputs measured in the experiments (real or virtual). • Solver of the direct problem (in most cases of high computation cost). • Optimization procedure of objective function. The objective function is defined as a distance between measured and calculated model outputs in the selected space norm or it can be Pareto set. Either gradient and non-gradient or bio-inspired optimization algorithms are applied to determine the minimum of the objective function.  Figure 1: Inverse analysis flow chart. 2.2 Sensitivity analysis Sensitivity analysis allows to assess the accuracy of the model of the analyzed system or process, determine the parameters which contribute the most to the output variability, indicate the parameters which are insignificant and may be eliminated from the model, evaluate these parameters which interact with each other, determine the input parameters region for subsequent calibration space [7],[8]. The steps of the sensitivity analysis are as following: • Sensitivity measure. The measure expresses the model solution (model output) changes to the model parameter variation.  • Selection of the parameter domain points. Design of experiment techniques are commonly used to select the lower number of points guaranteed searching whole the domain. • Method of sensitivities calculation. The sensitivities are estimated by global indices or by local ones. The information obtained from sensitivity analysis is applied to the inverse method: • To verify if the objective function is well defined – it means if it is possible to estimate the 1362D. SZELIGA  4parameters, which are looked for, based on the information included in the objective function. In case of no sensitivity or low sensitivity of the objective function to the parameter changes, the parameter identification cannot be performed and the objective function hast to be transformed to another form including verification of the model output space norm. • As the preliminary step – to select the starting point/the first region of interest or the first population for optimization algorithm. • Optimization process – to construct the hybrid algorithms (e.g. the combination of a genetic algorithm to select local minima and a gradient method to explore those minima) or modified algorithms (e.g. the particle swarm procedure enriched with the local sensitivities information [9]) to increase the procedure efficiency.  In this work the model output was defined as the Euclidean distance to exact logging curve. To points selection Latin hypercube  sampling (LHS) was used. The sensitivities were defined as the first order local sensitivities estimated using partial derivatives.  2.3 2D DC borehole resistivity measurements problem Computational domain. The problem geometry was described as 2D problem of plane coordinates (x,y). The following materials were used (figure 2a):  • borehole: a subdomain 0Ω  of width 10 cm ( ){ }0 , : 0cm 10cmx y xΩ = ≤ ≤  with resistivity 0 0.1 mρ = Ω⋅ , • upper and lower formations (no. 1 and 4): a subdomains 1 5,Ω Ω  defined by ( ){ }1 , :10cm , 3mx y x yΩ = < ≤ , ( ){ }5 , :10cm , 2mx y x yΩ = < < −  with resistivity 1,4 1000 mρ = Ω⋅ , • formation no. 2: a subdomain 2Ω  defined by ( ){ }1 , :10cm , 2m 3mx y x yΩ = < ≤ <  with resistivity 2 5 mρ = Ω⋅ , • formation no. 3: a subdomain 3Ω  defined by ( ){ }3 , :10cm , 0m 2mx y x yΩ = < ≤ <  with resistivity 3ρ , • formation no. 4: a subdomain 4Ω  defined by ( ){ }4 , :10cm , 2m 0mx y x yΩ = < − ≤ <  with resistivity 4 1 mρ = Ω . Variational problem formulation. Find Vu∈  the electrostatic scalar potential such that: ( ) ( )( )( )2121,,Ni i ii ib u v l v v Vu vb u v dxx xJl v v dS g v dSxσ=Ω=Ω Γ= ∀ ∈∂ ∂=∂ ∂∂= +∂∑∫∑∫ ∫ (6) where 1363D. SZELIGA  5( ) ( )2 22 : : tr 0 on DV v L v v dx vΩ⎧ ⎫= ∈ Ω + ∇ < ∞ = Γ⎨ ⎬⎩ ⎭∫  (7) and J denotes a prescribed, impressed current source, σ is the conductivity, and the electrostatic scalar potential u is related to the electric field E by E u= −∇ . More information of DC borehole resistivity measurements problem is presented in [10]. The direct problem formulated by equations (6) is solved using automatic hp-Finite Element Method software (for detail description see [11],[12]). An example of the hp adaptive computations is presented in figure 2b.   Figure 2: a) The 2D geometry of 3D DC borehole resistivity measurements problem. The rock formation is composed of five different layers of various resistivities, b) Automatic hp-adaptive solver – a solution example. 3 CALCULATIONS Inverse calculations. The objective function in the inverse method for DC borehole logging curve measurements was defined as the Euclidean distance between measured and calculated values: ( ) ( )211 mN i imi iu uN u=⎛ ⎞−Φ = ⎜ ⎟⎝ ⎠∑aa  (8) where u is the potential of electrical field, a – vector of identified parameters, N – number of measured points along the logging curve, m index – measured value.  The measurement values of electrical field potential um were generated for the problem ρ ?5 Ω⋅m1000 Ω⋅m1 Ω⋅m1000 Ω⋅myxborehole, 0.1 Ω⋅m2m1m2mΩ1Ω5Ω2Ω3Ω4Ω010cm1st step of the iterationInitial mesh(coarse mesh) Fine mesh Optimal meshp=1P=2P=3p=4p=5p=6p=7p=82nd step of the iteration1364D. SZELIGA  6described in the chapter 2.3 and the resistivity ρ3 = 200 Ω·m for Ω3 layer and next randomly perturbated. In general case vector a is of the form: { }, ,j jhρ α=a  (9) where ρj, hj – resistivity and height of the jth layer, respectively, α – deviation angel of the well from perpendicular,  j = 1…nf, nf – the number of formation layers. As yet inverse calculations were performed using hierarchical genetic searching [13],[14] with respect to some resistivities and the angel. Since the results of estimation were sufficient, the computations have been very time-consuming. The idea was to performed sensitivity analysis to investigate the parameters domain and to develop more efficient searching algorithm.  Calculation time depends on the assumed accuracy of hp-FEM solver. The accuracy of hp-FEM modeling is defined as the difference between coarse-grid and fine-grid solution in the quantity of interest. To increase the efficiency of the calculations some simulations of the process were performed with various accuracies. The calculations were carried out with the computer of two 3 GHz processors and 8 GB RAM, the times are shown in table 1. The accuracy of 10-2 was taken for further investigations.  Table 1: Relation between hp-FEM accuracy and the execution time of the solution. hp-FEM accuracy Execution time 10-5 5h 21min 40s 10-4 3h 30min 10-3 1h 4min 15s 10-2 20min 24s 10-1 9min 40s 1 3min 52s 10 2min 26s 100 2min 17s 1000 2min 19s  Sensitivity analysis. The sensitivity measure was expressed as derivatives of the objective function (8) and estimated though Taylor series expansion: ( ) ( ) ( )2 121nfii iaa+=∂ΦΦ +∆ = Φ + ∆ +Ο∂∑a a a a  (10) where ( )2Ο a  residue is neglected.  As the result the first-order local sensitivities are obtained { } { }/i is a= = ∂Φ ∂S . Current investigations were focus on the resistivity coefficient ρ and the sensitivity analysis was performed with respect to that parameter: / ρ∂Φ ∂ .  The domain for the resistivity was specified as [ ]0.1,1000ρ ∈  Ω·m. To provide uniform but random covering of the interval, points were generated by LHS with equal probability for 1365D. SZELIGA  7each subinterval. To evaluate numerical computations of / ρ∂Φ ∂  the value of ∆ρ was assumed as 5% and 10% of ρ resistivity. In figure 3 the logging curves calculated for various values of resistivity are presented. The changes of the goal function are shown in figure 4a and the results of sensitivities estimations – in figure 4b.  Figure 3: Logging curves of DC borehole resistivity problem: a) all formers layers, b) focus of one of the layer with curves obtained for various resistivity parameters.  a)  b)  Figure 4: a) Objective function computed for various resistivities, b) Sensitivity of the objective function with 0 0.2 0.4 0.6 0.8u, V-4048y, mρ1 = 100 Ω mρ2 = 5 Ω mρ3 = 200 Ω mρ4 = 1 Ω mρ5 = 1000 Ω m0 0.05 0.1 0.15 0.2 0.25u, V01234y, mResistivityρ = 1034 Ω mρ = 200 Ω mρ = 50.6 Ω mρ = 18.7 Ω mρ = 1.54 Ω mρ = 0.18 Ω m.a) b)0 200 400 600 800 1000 ρ, Ω m00.10.20.30.40.5Φ, −0 200 400 600 800 1000 ρ, Ω m-200204060ρ∂Φ∂1366D. SZELIGA  8respect to the resistivity parameter. 4 DISCUSSION Computations of the inverse calculations objective function sensitivities were performed for the resistivity screening the parameter domain. For resistivity values higher than 50 Ω·m the electrical field potential changes are not significant (see figure 3b). Similar behavior is observed for the objective function (figure 4a): the values are high for small resisitivties, close to zero for exact solution (200 Ω·m) and next they a little increase but are still low. The sensitivities are close to zero all over the interval except narrow interval of the exact solution (figure 4b). Such distribution of the objective function and the sensitivities effect the optimization procedure hard to performed using conventional optimization algorithms.  The investigations were performed for the resistivity of the one layer. In real systems the identification consists of several resistivities and the heights of the layers (vector a, equation (9)), consequently the objective function (8) is multimodal function of high sensitivities with respect to the resistivities close to the problem solution, remaining area is of low sensitivities.  The results of current investigations are guidelines for two optimization strategies: • Hierarchical genetic searching (or another bio-inspired optimization algorithm) to select local minima and next application of the gradient method to explore the local minima, • Application of design of experiment method (LHS or screening design) to select a set of points, for which objective function is computed with hp-FEM model and next, based on the those accurate results, generating fast metamodel using one of the approximation method (e.g. response surface algorithm) or neural network, and performing inverse calculations with investigated, fast model. Both the strategies are expected to increase the efficiency of the inverse calculations in terms of computation time. 5 CONCLUSIONS In that work application of the sensitivity analysis to the hard inverse problem as the preliminary step of the computations was proposed. As the example inverse problem of DC borehole resistivity measurements was formulated. The sensitivities were defined locally and expressed using Taylor series expansion. The set of the points that sensitivities were estimated for, was generated with design of experiment algorithm. The guidelines for the inverse optimization procedure were developed. The accuracy of the hp-FEM solution of the direct model in relation to the inverse problem objective function was analyzed as well. The obtained results of calculations, as the first attempt, are the basis to develop efficient optimization procedures with sensitivity analysis application to solve hard inverse problem: hybrid methods or algorithms with metamodels. Another way to decrease the number of the solver evaluations is to modify the objective function by including the information of the sensitivities. All above aspects will have been investigated. Acknowledgements. The work on sensitivity analysis has been supported by Polish MNiSW, project no. N N508 629 740. The hp-Finite Element Method calculations have been supported by Polish MNiSW, project no. NN519 447 739. 1367D. SZELIGA  9REFERENCES [1] Kirsch A., An Intoduction to the Mathematical Theory of Inverse Problems. Springer-Verlag, New York (1996). [2] Engl H.W., Hanke M., Neubauer A., Regularization of Inverse Problems. Kluwer Academic Publishers (1996). [3] Szeliga, D., Pietrzyk, M., Identification of rheological models and boundary conditions in metal forming, International Journal of Materials & Product Technology (2010), 39: 388–405. [4] Gawąd, J., Kuziak, R., Madej, Ł., Szeliga, D., Pietrzyk, M., Identification of rheological parameters on the basis of various types of compression and tension tests, Steel Research International (2005), 76: 131–137. [5] Szeliga, D., Gawąd, J., Pietrzyk, M., Inverse analysis for identification of rheological and friction models in metal forming, Computer Methods in Applied Mechanics and Engineering (2006), 195: 6778–6798. [6] Hadamard, J., Lectures on the Cauchy Problem in Linear Partial Differential Equations. Yale University Press, New Haven (1923). [7] Kleiber, M., Antunez, H., Hien, T.D., Kowalczyk, P., Parameter Sensitivity in Nonlinear Mechanics. Wiley (1997). [8] Saltelli, A., Chan, K., Scott, E.M., Sensitivity Analysis. Wiley (2000). [9] Sztangret, Ł., Szeliga, D., Kusiak, J., Sensitivity analysis as a supporting procedure in optimization of metallurgical processes, Hutnik Wiadomości Hutnicze (2010) 76(12): 721–725 (in Polish) [10] Pardo, D., and Torres-Verdin, C., Paszyński, M., Simulation of 3D DC borehole resistivity measurements with a goal-oriented hp-Finite Element Method, Part II: Through Casing Resistivity Instruments, Computational Geosciences (2008) 12(1): 83-89. [11] Demkowicz, L., Computing with hp-Adaptive Finite Elements. Vol. I, Chapman & Hall /CRC (2006). [12] Demkowicz, L., Kurtz, J., Pardo, D., Paszyński, M., Rachowicz, W., Zdunek, A., Computing with hp-Adaptive Finite Elements. Vol. II, Chapman & Hall/CRC (2007). [13] Barabasz, B., Migórski, S., Schaefer, R., Paszyński, M., Multi-deme, twin adaptive strategy hp-HGS, Inverse Problems in Science and Engineering (2011), 19: 3–16. [14] Schaefer, R., and Kołodziej, J., Genetic search reinforced by the population hierarchy, Foundation of Genetic Algorithm 7, De Jong, K., A., Poli, R., Rowe, J., E., Ed., Morgan Kaufman Publisher (2003): 383-399.  

Application of sensitivity analysis – preliminary step of the process parameters estimation

https://upcommons.upc.edu/bitstream/2117/184285/1/COMPLAS-2011_128-Application%20of%20sensivity%20analysis.pdf

Application of sensitivity analysis – preliminary step of the process parameters estimation

Abstract

Similar works

Full text

Available Versions

UPCommons. Portal del coneixement obert de la UPC