Presentation Outline

Source Characterization of Atmospheric Releases using Stochastic Search and Regularized Gradient Optimization Bhagirath Addepalli1, Christopher Sikorski2, Eric R. Pardyjak1, and Michael Zhdanov3 1Department Of Mechanical Engineering, University of Utah 2School of Computing, University of Utah 3Department Of Geology and Geophysics, University of Utah

Presentation Outline • Motivation • Problem Definition • Solution Methodology • Results & Discussion • Conclusions • Future Work

Motivation • The BioWatch Program • Department of Homeland (DHS) Security Program • Established in 2003 • Operational in 31 cities nationwide • Air samples tested daily for 6 biological agents of top concern • Spatial coordinates of the sensors known • If any of the sensors record a hit, can we track back to the source? Mock deployment of sensors in California

Problem Definition • The Forward Problem: A – Forward modeling operator (Gaussian plume model – GPM) m – Model parameters (Source parameters) d – Data (Concentration measurements at receptors) GPM for steady, continuous point releases under uniform wind conditions is: • Subscript ‘S’ represents source parameters • Subscript ‘R’ represents receptor parameters • Origin of the coordinate system at the source. The x, y, and z axes in the along-wind, cross-wind and vertical directions • σyandσz - Plume spread parameters – Account for turbulent diffusion of the plume

Problem Definition contd… • The GPM has eight model parameters • Five out of the eight model parameters retrieved in the present work. Hence the inverse problem is five dimensional (5D) • Values of C1 and C2 are terrain and problem dependent. C1 = 0.12 and C2 = 0.10 chosen for the data against which the proposed approach is validated

Problem Definition contd… • The Inverse Problem: • Model parameters (m) from the Copenhagen Tracer Experiments (TCTE) retrieved using the proposed approach θS=308.57° q S/ uS =0.64 g/m (0 m,0 m,115 m) Residential site

Inversion - Background • Operator A for the current problem is nonlinear. Hence, the problem is an “error minimization problem” (or) optimization problem • Derivative-based optimization adopted to take advantage of the simple analytical nature of GPM (analytical expressions for Frechet and Hessian can be pre-computed) • Solution to the inverse problem written as solution to the parametric functional: Regularization parameter Misfit functional Stabilizing functional • Gradient methods require the misfit functional to be convex, continuous and differentiable (C-C-D) • Misfit functional for the current problem is locally convex around the true source parameters • Stochastic search done to provide the gradient scheme a good initial iterate

Hypothesis • For the atmospheric inverse–source problem, the misfit functional looks like: Stochastic search algorithm Gradient Alg. Multiple Local Minima • The misfit functional has multiple critical points and a region in which it is convex shaped with a global minimum • CONCLUSION:To solve the source-inversion problem, a technique that comprises of stochastic search and gradient methods should be employed for fast reconstruction Global Minimum

Hypothesis contd… Illustration: region of convexity with global minimum To show the existence of the region of convexity and the global minimum, the TCTE problem was solved in two-dimensions by discretizing the domain and performing exhaustive grid search. The misfit functional based on L 2-norm was used ( )

Hypothesis contd… Illustration: region of convexity with global minimum Note: The extent (size) of the convex shaped region is a function of the domain size chosen

Solution Methodology L2-norm misfit functional Quadratic line-search Adaptive regularization Tikhonov stabilizer Solution Procedure log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search Gradient descent (Newton's method) Maximum sensors satisfied Final solution Trade-off between residual and sensors satisfied

Solution Methodology contd… Solution Procedure log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search • Done to account for zero measurements • Atmospheric inverse-source problems comprise of sparse number of measurements in general, and very few non-zero measurements in particular

Solution Methodology contd… Solution Procedure log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search • Not satisfying value at a sensor implies predicting a non-zero concentration at the sensor which in reality measured zero concentration (over-prediction) • Not satisfying value at a sensor implies predicting a zero concentration at the sensor which in reality measured non-zero concentration (under-prediction) • Limit of estimation (LOE) – the minimum concentration value a sensor can measure • LOE for Copenhagen tracer experiments LOE ≈ 9 × 10-9 g / m3 • ‘ε’ value was set at 10-16

Solution Methodology contd… Solution Procedure log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search • Rigorous strategy:“satisfy” the concentration measurements at all sensor locations. Cannot be implemented for all real-life dispersion problems • For the Copenhagen experiments, maximum number of sensors “satisfied” (Ns)= 131, out of the available 137 measurements (N). 1.4 × 109 MC simulations run for this exercise. Expected number of MC (E(MC)) points required for Ns= 131 was 9,790,210 (≈ 107) • Semi-rigorous strategy: “satisfy” “some” percentage of total number of measurements. Minimum value of NS for correct initial iterates = 128 (based on 1.4 × 109 MC simulations) • Rigorous semi-rigorous strategy: “satisfy” all non-zero measurements. For zero measurements, satisfy a percentage in the range [90%, 100%] – over-prediction

Solution Methodology contd… Solution Procedure log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search • Deterministic QMC point-sets are advantageous for atmospheric inverse problems because of their superior “space-filling” nature • Halton, Hammersley, NiederreiterXing, and SpecialNiederreiter point-sets taken from SamplePack (Kollig and Keller). Faure permutations applied to obtain the scrambled sequences • Sobol generator taken from Matlab. The MatousekAffineOwen scrambling procedure applied to get the scrambled form

Solution Methodology contd… Solution Procedure log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search • Stochastic (QMC) search performed in five dimensions (5D) • QMC points in the 5D hypercube [0, 1)5 linearly mapped to model parameter space. The domain of the model parameter space was:

Solution Methodology contd… L2-norm misfit functional Quadratic line-search Adaptive regularization Tikhonov stabilizer log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search Gradient descent (Newton's method) Solution Procedure • Newton’s method implemented with the initial iterate provided by the stochastic search stage • Value of q was set at 0.7 • Newton’s method took less than 100 iterations to converge to the global minimum for all tested cases

Solution Methodology contd… L2-norm misfit functional Quadratic line-search Adaptive regularization Tikhonov stabilizer Solution Procedure log10 (dobs / dpr) based misfit functional Deterministic QMC point-sets Semi-rigorous strategy Bounds for dobs / dpr Stochastic (QMC) search Gradient descent (Newton's method) Maximum sensors satisfied Final solution Trade-off between residual and sensors satisfied • ‘model-fit’ vs. ‘data-fit’ • In Newton’s method, should one go with predicted model parameters that minimize the residual or maximize the number of sensor measurements satisfied?

Results & Discussion • To get an estimate for the expected value of the number of MC points (E(MC)) (random samples) required to predict model parameters that satisfy a given number of sensor measurements, 1.4 × 109 MC simulations were run • Done for the semi-rigorous and rigorous semi-rigorous strategies N = 137; NNZ = 34; NZ = 103 • Semi-rigorous strategy: • min(NS ) = 128 – for initial iterates to be in the convex shaped region with probability 1 (apriori study) • E(MC) = 40,656 • Rigorous semi-rigorous strategy: • NS-NZ = 34; NS-Z = 94; NS = 128 • E(MC) = 439,147

Results & Discussion contd… • QMC point-sets used in the stochastic search stage • Number of points required for the semi-rigorous strategy (SRG) (min(NS)= 128) sensor measurements compared against E(MC) = 40,656 • Scrambled version of Halton point-set requires minimum number of points to satisfy NS = 128 • # of scrambled Halton points required = 797

Results & Discussion contd… • Number of points required for the rigorous semi-rigorous strategy (RSRG) (NS-NZ = 34, NS-Z = 94, NS = 128) sensor measurements compared against E(MC) = 439,147 • Original version of Halton point-set requires minimum number of points • # of original Halton points required = 3749

Results & Discussion contd… • Number of points required for the semi-rigorous strategy (RSRG) (min(NS) = 128) and NS = 129 sensor measurements compared against E(MC) = 92,075 • Scrambled version of SpecialNiederreiter point-set requires minimum number of points • # of scrambled SpecialNiederreiter points required = 15,505

Results & Discussion contd… • Number of points required for the rigorous semi-rigorous strategy (RSRG) (NS-NZ = 34, NS-Z = 95, NS = 129) sensor measurements compared against E(MC) = 503,778 • Original version of Sobol point-set requires minimum number of points • # of original Sobol points required = 47,213

Results & Discussion contd… • Newton’s method takes about 70 iterations to converge to the final solution when the semi-rigorous strategy (SRG) is employed - NS = 128 Max. sensors satisfied = 130 Trade-off between residual & sensors satisfied

Results & Discussion contd… • Newton’s method takes about 10 iterations to converge to the final solution when the rigorous semi-rigorous strategy (RSRG) is employed - NS = 128; NS-NZ = 34

Conclusions • Approach comprising stochastic search and regularized gradient optimization developed to solve the inverse-source problem • A new misfit functional and two strategies proposed for solving the inverse problem and identifying the initial iterate • Based on the tests conducted, the original versions of Halton and Sobol, and scrambled versions of Halton and Hammersley consistently perform better than the Mersenne-Twister generator

Future Work • Simulated Annealing (SA) to be employed to retrieve all the model parameters (8D problem) • Hybrid SA comprising homogeneous and inhomogeneous components will be implemented • QMC point-sets will be used in the homogenous stage • Will be run on GPU’s to obtain regions that contain the source parameters rather than a single model that fits the data

Presentation Outline