Surrogate models and optimal design of experiments for chemical - - PowerPoint PPT Presentation

CLEAN COMBUSTION RESEARCH CENTER

Surrogate models and optimal design of experiments for chemical kinetics applications

F . Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone | January 7, 2015

King Abdullah University of Science and Technology | Reactive Flow Modeling Laboratory

www.kaust.edu.sa

Fossil fuels sustain worldwide energy demand through 2030

& ( / + '% '* '. '.,& '+'& '+*& '++& %&(& 0121345617 89:614; <=>;? @47 AB6 C?46 "B66B?2$L?1

S$T2:69>17$5B?M9167 S Source: BP Energy Outlook 2030 (2011)

Energy growth is mainly related to GDP Largest growth in Non-OECD countries Pollutants (CO2 emissions, particulates, etc.) act as a constraint

The role of combustion

Combustion technology enables more efficient and cleaner processes: power, transportation, industrial

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 2/28

The complex network of reactions...

Hydrocarbon fuels (CnHmOp) and oxygen react to produce H2O and CO2 (and traces of pollutants...) with intense energy release (“heat”) This process occurs through thousands of intermediate species and reactions (e.g., H + O2 −

− ⇀ ↽ − − OH + O)

Overall rate of reaction controls: (a) heat release, pressure rise, ...; (b) stability and robustness of combustion process; (c) pollutant emissions

Chemical kinetic models

The formulation of a (complex) network of reactions involving many chemical species is key to the accurate and faithful description of combustion processes

⇒ Reaction are characterized by reaction rate parameters,

which are the model parameters we are concerned with in this talk

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 3/28

A single reaction and its description

Consider the chain branching reaction

H + O2 −

− ⇀ ↽ − − OH + O

The rate at which this reaction occurs (mol cm−3 s−1) is given by the following expression q = kfCHCO2 − kbCOHCO (1) where kf,b are the forward and backward rate “constants” and Ci the species concentration (mol cm−3). The “Arrhenius form” is customarily used to reflect T dependence kf = kf(T) = AT b exp(−E/RT) (2) kb = Keq(T)kf (3) A (pre-exponential factor), b (temperature exponent), and E (activation energy) are the model parameters θ = (A, b, E).

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 4/28

A fundamental tool in combustion: shock tube

driver (inert gas) driven section (test mixture) end wall diaphragm

• ptical

detector pressure sensors reflected shock shock T0, p0

Adapted from Petersen Combust Sci Tech 2009.

Allows to bring a gas to p0 and T0 instantaneously via a (reflected) shock Time evolution of spatially homogeneous mixture is modeled as system of ODEs Ideal for kinetics studies when complemented with advanced diagnostics:

τign, XOH(t), XCH4(t), etc.

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 5/28

Seeking to measure reaction rate parameters

• Goal. Apply a framework of optimal experimental design to the characterization
• f elementary reactions of interest to Prof. Farooq

H + O2 −

− ⇀ ↽ − − O + OH

Methylfuran + OH −

− → products

Key features Bayesian inference approach A measure of information gain A fast estimate of the expected information gain Polynomial Chaos expansion as surrogates for the model Sparse adaptive pseudo-spectral projection to build surrogates

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 6/28

Measuring the rate of H + O2 −

⇀ ↽ − O + OH (I)

Data gathering A shock tube is filled with H2, O2, Ar (inert diluent) in prescribed proportions The mixture is shock-heated at p0 and T0 and water H2O is formed The concentration of H2O is measured by time-resolved tunable diode laser absorption The ODE system is solved adjusting the rate constant until a good fit is found

0.0 0.2 0.4 0.6 0.8 1.0 500 1000 1500 2000 2500 Best fit k1 Best fit k1 x 1.1 Best fit k1 x 0.9 Experiment

(a)

Water conc. [ppm] Time [ms]

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 7/28

Measuring the rate of H + O2 −

⇀ ↽ − O + OH (II)

Selecting the composition Typically tests are run for very diluted fuel/rich mixtures: e.g., O2 at 0.1%, H2 at 0.9% and balancing Ar For these mixture compositions, the peak rate of water formation is very sensitive to H + O2 −

− ⇀ ↽ − − O + OH $

1 5 10 15 20 0.5 1 Reaction ID

• Sens. Index

1 5 10 11 20 Reaction ID Sensitivity Index 1 0.5 H + O2 = O + OH

• Note. The selection of “optimal mixture composition” is a result of years of

experience and considers issues of signal-to-noise ratio, detectability thresholds, and time resolution

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 8/28

Measuring the rate of H + O2 −

⇀ ↽ − O + OH (III)

Obtaining A, b, and E Tests are performed for various values

• f T0

For each test, the “best fit” k is found by solving the ODE system with perturbed values of k The pairs (k, T0)i are “best fit” to the Arrhenius expression

K 0.3 0.5 0.7 0.9 10

11

10

12

10

13

2x10

13

3333 K 2000 K

k1 [cm

3mol

• 1sec
• 1]

1000/T [1/K] Hong et al. (2010) Masten et al. (1990) Pirraglia et al. (1989) Hong et al. (2010)

1429 K 1111 K 4x10

10

Quoting from Hong et al.: “The rate is given as k = (1.04 ± 0.03) × 1014exp[(7705 ± 40)/T] cm3 mol−1 s−1 over the temperature range 1100 to 3370 K based on forward propagation of uncertainties

• n other reactions, accuracy of experimental diagnostics, knowledge of T0, etc.”
• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 9/28

Bayesian experimental design

We begin with a model for the experiment y = G(θ, ζ) + ǫ (4) y: observable (i.e, maximum rate of increase of [H2O])

θ: model parameters (Arrhenius coefficients of target reaction, A, b, and E) ζ: experimental design variables (e.g., temperature and mixture

composition) G: combustion simulation model (approximated by surrogate).

ǫ: measurement noise, assumed to be Gaussian.

The evaluation of G requires the solution of an initial value problem featuring a system of ODEs of relative large size

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 10/28

Information gain metric

We adopt the Kullback-Leibler divergence as information gain DKL =

• Θ

log p(θ|y) p(θ) p(θ|y)dθ. (5) with the expected information gain, I = E[DKL] =

• Y
• Θ

log p(θ|y) p(θ) p(θ|y)dθp(y)dy. (6) where p(θ): prior p(θ|y): Bayesian posterior

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 11/28

Computing the expectation of DKL

The computation of I = E[DKL] is challenging as it must be computationally inexpensive if seeking to maximize I as a function of experimental conditions ζ. Double loop Monte Carlo Laplace approximation to DKL We use a Laplace approximation (Long et al. 2013, 2014) Analytical integration against the posterior measure The calculation of I is reduced to a single loop Requires derivatives of the model (Hessian matrix) Because it’s fast, it enables searches for the optimized experimental designs

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 12/28

Surrogates

Fast model surrogates for the estimations of DKL and I are required and truncated Polynomial Chaos expansion (PCE) are used to construct such surrogates G(ξ) ≈

N

• k=0

ckΨk(ξ) (7)

ξ = (ξ1, ..., ξα) : a vector of germs which maps model parameters θ and

design variables ζ. G(ξ): Quantities of interest (QoI).

{Ψk}: a set of multi-dimensional Legendre polynomials.

ck = G, Ψk/Ψk, Ψk obtained via sparse adaptive pseudo-spectral projection (Winokur et al. 2013) The surrogates provide derivatives needed for Laplace approximation

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 13/28

Estimation of H + O2 −

⇀ ↽ − O + OH

The goal is the estimation of (1) A and (2) E for the forward rate constant (assume b = 0) The design parameters are (1) H2 concentration; and (2) temperature T0

⇒ [H2]0 ∈ [0.5 − 5]% ⇒ T0 ∈ [1100 − 1500] K

Pressure is constant at p0 = 2 atm and [O2] = 0.1%. The observable is max{d[H2O]/dt} during the mixture evolution Mapping for each random model parameters is xi = xi,0 + xi,0diξi

(i = 1, 2) ξi ∼ U[−1, 1]

(8) The nominal values A0 and E0 are from Hong et al. (2011) and complemented by dA = 0.5 and dE = 0.1, respectively (giving a change in the rate constant from −44% to +80%)

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 14/28

Pseudo-spectral Projection: Sensitivities and Surrogates

Global sensitivity indices are computed with 20 reactions with ±50% uniform uncertainties on ki (i = 1, . . . , 20) Reaction 1 (H + O2 −

− ⇀ ↽ − − O + OH): sensitivity index > 0.98

A/A0 ∼ U[0.5, 1.5] & Ea/Ea,0 ∼ U[0.9, 1.1]

Mean and variance

0.5 1 0.05 0.1 0.15 0.2 time(ms) [H2O](%) µ µ±σ

Surrogate

0.5 1 1.5 5 10 15 A/A0 max slope 0.9 0.95 1 1.05 1.1 5 10 15 Ea/Ea,0

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 15/28

Expected information gain for 1 exp

Let us focus on a single experiment scenario first

1 G-L quadrature point 5, 13 G-L quadrature point

I([H2]0, T0) computed with Laplace approximation converges with ≥ 5 quadrature points I depends strongly on T0 (higher for higher T0) I has a weak dependence on [H2]0, except for very low H2 concentrations

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 16/28

Non-optimal and optimal 1 exp design

Based on the I surface, we consider two single experiments – A “non-optimal”: low H2 concentration and low temperature – The “optimal”: high H2 concentration and high temperature

1100 1200 1300 1400 1500 2 4 6 8 10 12 T(K) growth rate(ms−1) 0.5 1 1.5 0.9 0.95 1 1.05 1.1 A/A0 Ea/Ea,0 DKL=2.3 DKL=1.0

The posterior (sampled with a Markov Chain Monte Carlo) distribution is narrower for the optimal design and the corresponding DKL is higher

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 17/28

Some considerations

Disclaimer: the “one experiment” scenario leads to an “underdetermined” problem for the two model parameters A and E Posterior arranged around two different k = const curves on the A–E plane (recall the two temperatures are different giving different k’s) The larger scatter for the lower T case most likely reflects the smaller signal to noise ratio An explanation: “optimized” design really seeks to maximize the signal to noise ratio (σǫ = const) by increasing the temperature, thereby leading to a greater rate of formation of water (observable) The I surface suggests best working with rich mixtures, in agreement with “expert advice” (rich mixtures for [H2]0 > 0.2%)

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 18/28

Expected information gain for 2 exp

Now consider a two experiments scenario whereby the concentration of H2 is fixed and two initial temperatures are chosen. Control parameters: T0,1 and T0,2 [H2]0 = 5.0% fixed IKL,max at T0,1 = 1226 K and T0,2 = 1500 K The surface is symmetric as expected Maximum is along the T = 1500 K plane Minimum is along the T1 = T2 surface

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 19/28

Executing 20 experiments

Now execute 20 runs, according to the following principles Scenario 1 (“optimal”): 10 runs at 1500 K and 10 runs at 1266 K Scenario 2 (“uniform”): 20 runs at equi-spaced temperatures

1100 1200 1300 1400 1500 2 4 6 8 10 12 T(K) growth rate(ms−1) 0.5 1 1.5 0.9 0.95 1 1.05 1.1 A/A0 Ea/Ea,0 DKL=4.4 DKL=3.9

1.1 0.9 0.95 1 1.05 1.1 10 20 30 40 50 k/k0 prior design1 design2 1100K 1500K

The difference in DKL and posterior between the two sets is not large The “optimal” set provides a higher DKL and a narrower posterior and may be more convenient experimentally The conclusions reflect the model k = A exp(−E/RT) strongly

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 20/28

Some considerations...

The uniformly distributed temperatures (scenario 2) is more or less what is done today As part of the Bayesian inference process, a posterior distribution for A and E is obtained The posterior shows strong correlation between A and E and may be used for forward propagation UQ studies Contrast posterior with the independent, uniform random variables used for the prior (or as typical input for forward propagation UQ studies) Move beyond k = (1.04±0.03) × 1014exp[(7705±40)/T] cm3 mol−1 s−1

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 21/28

Measuring the rate of Methylfuran + OH −

→ products

The shock tube is filled with Methylfuran (MF), tert-Butyl hydroperoxide (TBHP) and Ar inert diluent TBHP serves as a OH precursor via thermal decomposition with a well-known rate OH radical concentration is measured via absorption spectroscopy

tert-Butyl hydroperoxide

20 40 60 80 100 5 10

t (µs) [OH] (ppm)

T=900K, p=1.5atm, [MF]=200ppm, [TBHP]=20ppm

τd The concentration of OH increases first (due to TBHP decomposition) and decreases next due to consumption by the target reaction MF + OH −

− → products

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 22/28

Experimental Design, Methylfuran

Observable: decay time τd of [OH] Control parameters: T0 (900 – 1400 K) and [TBHP]0 (10 – 20 ppm) MF + OH ↔ furylCH2 + H2O k = AT be−E/RT = (AT b

ref)(T/Tref)be−E/RT = A′(T/Tref)be−E/RT

A′ = AT b

ref , b, and E

Tref = TlowTup = 1122 K

0.7 0.8 0.9 1 1.1 10

11

10

12

10

13

10

14

10

15

A′/A0

′ : 5±

1000/T k 0.7 0.8 0.9 1 1.1 10

11

10

12

10

13

10

14

10

15

b/b0: ±100% 1000/T 0.7 0.8 0.9 1 1.1 10

11

10

12

10

13

10

14

10

15

Ea/Ea,0: ±100% 1000/T

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 23/28

Pseudo-spectral Projection: Surrogates

Methylfuran oxidation: 306 species and 1483 reactions A′/A′

0 = 5ξ, ξ ∼ U[−1, +1],

b/b0 ∼ U[0, 2.0] Ea/Ea,0 ∼ U[0, 2.0]

1 20 40 60 80 A′/A′ taud,[OH](µs) 1 2 20 40 60 80 b/b0 1 2 20 40 60 80 Ea/Ea,0

E0 Emin Emax A′ A′

min

A′

max

A′ A′

min

A′

max

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 24/28

Experimental design, Methylfuran

Double Loop Monte Carlo Imax = 3.17 @ Tlow & [TBHP]low Low concentrations of TBHP and low temperatures are preferred

900 1000 1100 1200 1300 1400 10 12 14 16 18 20 T0(K) [TBHP]0(ppm) 2.2 2.4 2.6 2.8 3

“Underdetermined”: surface (close to a plane) with scatter around it Maximization of signal to noise (noise is constant) appears to guide maximization of I, since T0 ↓, τd ↑. Minimization of TBHP is consistent with expert advice, but no trade-off with respect to detectability of OH is considered

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 25/28

Conclusions I

Successfully built surrogates for experimental models in reaction kinetics studies Surrogates are a key ingredient for forward propagation, inference, and

• ptimal experimental design, etc.

Coupled surrogates with optimal experimental design based on maximization of the expected information gain Two applications from ongoing experimental campaigns: chain-branching reaction & furan oxidation

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 26/28

Conclusions II

Some thoughts... The noise model “casts a long shadow”. Is the optimal design seeking to maximize the signal to noise ratio, or expose the sensitivities of the

• bservables to the model parameters to be inferred? Or both?

The choice of the observable is key and tightly connected to our ability of formulating a realistic noise model with a “manageable” functional form The ranges of design variables are somewhat “narrow” as they come from “expert advice”. Are the I surfaces “flat” because we are limiting ourselves to “obvious” design choices? Need to incorporate more realistic “sensor models” with associated noise

• models. Would this bring a more “interesting” I surface with “physics-related”

trade-offs leading to local minima and maxima?

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 27/28

Thank you!

fabrizio.bisetti@kaust.edu.sa http://flow.kaust.edu.sa

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 28/28

Equations

A homogeneous mixture of ideal gases is described by z = (p, T, Y1, . . . , YM) T is temperature and p is pressure Y = (Y1, . . . , YM) are the mass fractions of species (M = 10 − 1000): e.g. YO2 = mO2/m = ρO2/ρ. Yields a system of ODEs of size M + 2

∂Yi ∂t = Wi ρ ωi

i = 1, . . . , M (9)

ωi = ωi(T, p, Y) =

R

• r

(ν′′

r,i − ν′ r,i)qr

(10) qr = qr(T, p, Y) = kr

f (T) M

• i

C

ν′

r,i

i

− kr

b(T) M

• i

C

ν′′

r,i

i

(11)

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 29/28

Expected information gain I

Double loop Monte Carlo (DLMC) I ≈

NI

• i=1

log p(yi|θi) p(yi) (12) p(yi) ≈

NJ

• j=1

p(yi|θj). (13)

Simple method to estimate the expected information gain I. Not fast enough to search for an optimal experiment.

Markov chain Monte Carlo (MCMC)

Efficient statistical inference technique by sampling the Bayesian posterior for a given experiment data y. Posterior samples for further analyses.

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 30/28

Expected information gain II

Estimate the information gain DKL DKL ≈

NMCMC

• i=1

log p(θi|y0) p(θi) , (14)

{θi} is an MCMC posterior sample for given experiment data y0.

Not fast enough to estimate the expected informatio gain, I in reasonable computation times.

Fast estimation of expected information gain based on the Laplace approximation

Enables searches for the optimized experimental designs. Requires derivatives of the model.

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 31/28

Expected information gain III

Sensitive to measurement errors, Σǫ. DKL = −1 2 log

• (2π)d|Σ|
• − d

2 − Σ∇∇h(ˆ

θ)

2

+ OP

• 1

M2

• ,

(15) where

Σ = H−1

and H =

M

• i=1

∇G⊤

i Σ−1 ǫ ∇Gi

(16)

• F. Bisetti, A. Farooq, D. Kim, O. Knio, Q. Long, R. Tempone – Surrogate models and optimal design of experiments for chemical kinetics applications

Jan 7, 2014 32/28