Local and global parametric analysis of reacting flows

Physica 20D (1986) 67-81 North-Holland. Amsterdam

LOCAL AND GLOBAL PARAMETRIC ANALYSIS OF REACTING FLOWS Herschel RABITZ Department

of Chemistry, Princeton University, Princeton, NJ 08544, USA

Received 10 December 1984 Revised manuscript received 29 January 1985

This paper presents a selected overview of current techniques available for analyzing parametric behavior of dynamical systems including reactive flows. The paper is based on a series of lectures presented at Cornell University, and the discussion is broadly divided into the reahns of local and global techniques. Local procedures are capable of addressing a rich variety of questions associated with the behavior of dynamical systems near a chosen operating point in parameter space. On the other hand, global procedures are specScally designed to map out solution behavior over wide ranges of parameter space for both statistical analysis as well as objective function mapping. Techniques for accomplishing all of these tasks are in various stages of development and comments on current research directions are mentioned when appropriate.

1. Introduction The development and application of numerical techniques for solving differential equations describing reactive flows continues to be an .active area of study. Computational considerations will surely remain to be an important aspect, especially for complex systems. The analysis of the solutions to such problems with regard to their parameter content has also received some attention over the years. However, with the increasing availability of sophisticated numerical codes for solving the differential equations, there is becoming increasing interest in exploring solution behavior as a function of the underlying model structure or parameters. Historically, the first area of concern in this regard was with respect to matters of stability analysis, and recently there has been a resurgence of interest in this topic. In addition, regardless of whether a dynamical system is stable or not, there is always the basic desire to know which portions of a model are responsible for particular observations. The subject of stability analysis and parametric control are, in fact, closely related, but for physical reasons they are often thought of as dis-

tinct. Stability analysis is usually associated with system response to initial condition variations. In the context of the present paper, the latter conditions are merely included amongst all the system parameters to produce a unified approach. The matters addressed above all fall within the context of sensitivity analysis [l]. There are two general reasons why these topics are of significance in the study of dynamical systems. First, many of the parameters, including initial conditions, are not precisely known in a given problem. Therefore, statistically related issues arise concerning how a given level of input uncertainty propagates to produce correlation or uncertainty in the output observables. Regardless of parametric uncertainty problems, there is always the overriding matter of the physical role played by each parameter in a system. Dealing with this latter question will ultimately be the topic of prime concern particularly as the input information becomes better defined in a given system. The study of these topics goes beyond merely solving the original dynamical equations, and additional cost can be involved. In accord with the varied nature of the questions above, a wide variety of sensitivity related tech-

0167-2789/86/$03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

68

H. Rabitr / Local and global parametric ana!vsis of reactingjlows

niques are continuing to be developed [2, 31. All of these procedures fall into either the class of local or global methods. The distinction between these two groupings can be understood by first recalling that the numerical solution to any model first requires specification of the system parameters at an operating point in parameter space . Therefore, a local analysis entails probing the solution behavior in the immediate vicinity of the chosen operating point. In contrast, global sensitivity analysis aims to map out solution behavior over wide regions of parameter space. Clearly, the latter topic will encompass the former one, however, the considerable disparity in the labor involved warrants treating them as separate subjects. Local techniques have been developed to a higher degree of sophistication and they go far beyond merely asking how a given observation is related to a particular parameter variation. For example, local techniques are capable of addressing the interrelationship amongst all the dependent and independent variables in a system. In addition, local techniques may be easily utilized to address the parametric origin of characteristic features in the observable output. Local sensitivity analysis techniques are under continuing development with regard to their numerical implementation and practical application. This ongoing effort has been currently reviewed rather extensively [2, 31, and the present paper will not give the topic special emphasis. The essence of the local procedures will be discussed in this paper with a focus on certain specialized topics that have not already been extensively pursued, especially in the areas of chemical dynamics and reactive flows. On the other hand, the topic of global sensitivity analysis appears to only be in its infancy and the ultimate directions of the field are not entirely clear at this time. The paper will, therefore, present an understanding of the subject such as it exists within the context of the physical problems addressed here. A special effort will be made to connect the global and local approaches whenever possible. Therefore, the following sections will be organized to treat local

methodology (section 2), global concepts (section 3), and finally, some general comment (section 4). Extensive numerical results will not be presented in this paper since they can be found in the current literature. In this regard, the paper also will not attempt to provide an exhaustive survey of the literature. Rather, selected references leading to the remaining literature will be cited where appropriate.

2. Local sensitivity analysis techniques In general, reacting flows will often be described as either initial value problems, boundary value problems, or mixed initial-boundary value problems. Although, in strict terms, steady boundary value problems do not show temporal evolution, they are still of interest in a dynamical sense. This comment follows since nominally steady systems may ultimately go through temporal excursions due to the presence of natural fluctuations. Sensitivity analysis in each of these classes of problems has its own particular characteristics. Nevertheless, there are certain similarities and the emphasis here will be from this perspective. This point of view is particularly relevant for two reasons. First, the practical methodology for calculating sensitivity coefficients is beginning to converge upon a common set of schemes showing similar structure regardless of the mathematical nature of the original physical problem. Secondly, once the raw sensitivity coefficients are obtained, their interpretation and further manipulation is entirely similar from one class of problem to another. 2.1. Numerical evaluation of sensitivity coefficients and densities Within the types of differential equations discussed above for physical modelling, there are two broad classes of sensitivity problems that arise corresponding to whether the parameters and their variations are constant or treated as distributed functions of space and/or time. In the former

H. Rabirz/ Local and global parametric analysis of reactingjows

case, local sensitivity analysis has largely been concerned with the partial derivative gradients 80i(r, t)/&xj of the observables Oi(r, t), i = 1, . . . , N at time t and position r with respect to the system parameters aj, j = 1,2,. . . , M. In the second category, when the parameters are distributed, it is necessary to extend this concept to the analogous functional derivatives SOi( r, t)/ Saj(r’, t’) where these latter densitiesare the basic response functions connecting incremental parametric deviations &Jr’, t’) to corresponding changes in the observables

(1) We shall focus here on Fr6chet derivatives, but analogous Gateaux derivatives could also be examined. Even when the system parameters and their variations are treated as constant, it is still informative to calculate functional variations since a more detailed understanding of the role of the system parameters can be achieved. Furthermore, in this latter circumstance, it may readily be shown that the familiar partial derivatives can be recovered by a simple integration. ‘Oitry

‘1 =

aa,

udr’dt’ SOi(r,t) Saj(r’,

2’) *

In the remainder of the paper on local sensitivity methodology, functional or partial derivatives, will be treated wherever appropriate. In order to be specific, we may consider a physical problem to be modelled by differential equations (it should be noted that equations of an integral or other type could also be similarly treated) of the following general form: L(O,a)

=o.

(3)

Eq. (3) represents a coupled system of equations which is typically non-linear in the output observables. Problems in reacting flows are naturally

69

of this type. In order to completely specify the physical model, the appropriate initial and/or boundary conditions would accompany eq. (3). All this information, including the nominal value of all the system parameters, is assumed to be known in order to proceed with the mathematical modelling process alone. In practice, of course, the model may be incomplete and secondly, the parameters may not be precisely known. This situation represents a standard realm of application for sensitivity analysis with the goal being understanding how uncertainties in the input map onto the observable output. Although local sensitivity gradients may be directly applied in statistical formulas for this purpose, ultimately the question becomes a global one when there are large uncertainties in the system parameters and discussion in this regard will be given in section 3. Sensitivity analysis at both the local and global levels represents a powerful set of tools regardless of the issue -of parameter uncertainty. For example, even if the parameters in a system are precisely known, their role in controlling the behavior of particular observables can be severely obscured by virtue of the convoluted “mixing” of the parameter contributions upon solving the system eqs. (3). For this reason, sensitivity analysis should best be thought of as a reliable means of quantitatively analyzing the parametric or physical content ensuing from a mathematical model of the particular physical process. A local sensitivity analysis begins by assuming that eq. (3) may be solved in a satisfactory fashion. An important point in this regard concerns the generally nonlinear nature of these latter equations, which therefore prescribes the need for iteration in the solution process. A limitless number of iteration schemes could be envisioned, but the Newton technique [4] is especially advantageous when also seeking sensitivity coefficients; this point will become clear below. A Newton iteration scheme based on eq. (3) would take on the following form J”( r, t)AO”( r, t) = -L (O”, a),

(4)

H. Rabitz/ Local and globalparameiric analysis of reacringjows

70

where A&” is the correction to the nth iterate 0” (r, t) with On+’ = 0” + AO”. The Jacobian J’(r, t) is a differential operator having elements Jjj(r, t) = 6’Lj( On( r, t), a)/dO~(r, t). Practical schemes based on eq. (4) may be generated including the use of adaptive (variable) space and time stepping. In order to obtain sensitivity information, we take the functional derivative of eq. (3) with respect to the system parameters +,

l)

wr,t’)0 = __&x(r’, WC t> t’) ’

8a( r’,

where an element of the sensitivity density matrix is understood to have the form

Wr, t)

afli(r, f)

&a( r’, t’) ii = &xj( r’t’)

’

Elements of the right-hand side of eq. (5) have a similar structure and the arguments of the operator L are given in terms of a point in space and time for notational convenience. The Jacobian differential operator J entering into eq. (5) is assumed to be evaluated at the exact or numerically precise solution to eq. (3). The evident parallel structure of eqs. (4) and (5) can be immediately exploited to obtain the sensitivity densities at virtually no extra cost over that of solving eq. (4). A practical algorithm would be attained by iterating eq. (4) until a satisfactory solution is achieved and then turning to eq. (5), at which point the Jacobian in the latter equation would be identified with that in eq. (4). The efficiency of this procedure for solving eq. (5) may be understood by realizing that the solution to eq. (4) effectively requires obtaining the inverse of the Jacobian operator. Generation of J and its inverse represents the bulk of the computational effort in solving eq. (4) and having J -’ operate on the additional inhomogeneity of eq. (5) incurs very little additional expense. The procedure presented here may be applied in various forms to a broad class of mathematical problems. In practice at this time the technique has only been implemented in a restricted fashion to

the initial value problems, one-dimensional boundary value problems, and strictly algebraic systems for sensitivity gradients rather than densities [5]. The procedure argued above made explicit use of joining the solution process for the original physical equations with that of the sensitivity equations. Assuming that sensitivities are to be desired along with solving the original model equations, any increased efficiency beyond that of the present prescribed algorithm must await the development of possibly better numerical procedures for the model eqs. (3). Finally, there is one point of caution concerning this procedure. In particular, due to the reliance on joining the solution and sensitivity procedures, a basic approximation is made. It is assumed that the sensitivities have essentially the same spatial and/or temporal variation as the system solution. This assumption is embodied in using the same discretized Jacobian in eqs. (4) and (5). From presently available results this central assumption seems to be generally valid, although serious exceptions could easily be envisioned.

2.2. Sensitivity

Green’s functions and their role

Although the inverse of the Jacobian operator was discussed above, in practice, such an inverse would likely not be explicitly calculated. Rather, a discretization of the operator would be achieved on a space and time grid followed by a decomposition or factorization (e.g., Gaussian elimination) to achieve a solution. Nevertheless, from a sensitivity point of view, it is physically important to calculate the inverse. In addition, the use of Green’s functions (equivalent to the inverse) have historically been utilized as an alternative means to solving eq. (5) and for the study of certain types of objective function sensitivity questions (cf., section 2.3). The physical significance of the Green’s functions may be understood by introducing an incremental source term into eq. (3) L(o,

a) = V(r,

t),

(6)

H. Rabitz/ L.ocal and globalparametric analysis of reactingflows

where it is understood that 6/&r, t) is an incremental flux associated with the observable Oi( r, t ). In the case of combustion, such observables would typically be species concentrations or the temperature. Functional differentiation of eq. (6) with respect to elements of the introduced flux vector will produce the following equation: J(r, t)G(r, t; r’,t’)

=lcY(r-r’)S(t-t’),

(7)

where Gii(r, t; r’t’) = 60i(r, r)/Q$(r’,

t’).

(8)

The connection of G and the inverse of J is evident from the structure of eq. (7). In practice, the numerical calculation of G would be achieved by merely solving eq. (5) once again with a unit inhomogeneity. For a purely temporal problem, one may also show that the Green’s function reduces to

Gij( I, t’) =

ao,(t)/ao&‘), 6ij, i 0,

t > f’, t = f,

(9)

t < t’.

The causality conditions in eq. (9) also apply to the more general space time Green’s functions in eq. (8). The physical significance of the system Green’s function derives from two sources. First, elements of the Green’s function matrices may be thought of as special sensitivity coefficients or densities. In particular, eq. (8) would be interpreted as the sensitivity of the ith observable at position r and time t with respect to a disturbance in the jth species at position r’ and time t’. Secondly, as stated above, the solution to the sensitivity equations may be explicitly expressed in terms of the system Green’s function. These two points will now be discussed separately. Except for problems where the parameters aj, j = l,l,. . . , M include laboratory control variables, the sensitivity coefficients or densities are not in themselves directly measurable quantities. With this exception aside

71

the utility of sensitivity coefficients or densities lies in their ability to indicate the role of the parameters in a model. In this sense, Green’s function elements as sensitivities are special since the elements of the Green’s function matrix may be measured in the laboratory. For example, one could introduce an incremental flux of a given species in a reacting flow and monitor the response of the entire pool of chemical species. In the limit that the disturbance is truly infinitesimal, the measurements would exactly correspond to the elements of the Green’s function matrix. By following the above line of reasoning, the magnitude of the Green’s function matrix elements may also be identified with the stability of the dynamical system. Indeed, the standard Liapunov stability numbers may be calculated by diagonalizing the Green’s function matrix. For example, in the purely temporal case of eq. (9) the eigenvalues A( t, t’) would be obtained by determining the similarity transformation U(t, t’) such that A&, t’) = [U-‘(t,

t’)G(t, t’)U(t, t’)lii

(10)

In transient systems, the evolutionary stability is often the question of concern rather than the asymptotic stability (e.g., a closed reactive system will always exhibit asymptotic stability due to thermodynamic considerations). Therefore, an interesting question is to compare t-’ In Xi(t, 0) with the eigenvalue ai of the Jacobian J(t). This comparison is significant since the Green’s function is determined by the Jacobian matrix through the formal relation G(t,t’)=

Texp[[:J(r)dr],

where T is the time ordering operator. Calculations have indicated that, in explosive systems, the eigenvalues to the Jacobian may underestimate the explosion time and overestimate the intensity of the explosion as compared to that prescribed by the Green’s function [a]. The elements of the transformation matrix U(t, t’) indicate which chemical

12

H. Rabilz/ Local and global parametric ana!ysis of reuctingjlows

species participate in’ the growth or decay of the system as time evolves. The second role for the Green’s function enters through its connection to the solution for the sensitivity eqs. (5)

Wr, t) &Y(r’, f’) = -

JJdt” dr” G( r, t; r”, ,“) (11)

In the case of purely temporal problems an efficient algorithm [7] has been written to calculate the sensitivity coefficients based on this integral expression, although the scheme embodied by exploiting the similarity of eqs. (4) and (5) is likely more efficient. An important point to again consider in this regard is that the latter procedure would lose its degree of efficiency if the sensitivity equations require discretization in a manner entirely different from that of the original model equations, and the alternative provided by calculating the Green’s function and utilizing eq. (11) could circumvent this problem when it occurs. 2.3. Objective functions and their sensitivities The observables discussed thusfar in a reacting chemical system are merely the direct output of the system differential equations. These would typically include the species concentrations, temperature, pressure, etc. However, in many realistic problems, other quantities may be of more direct interest. In this regard, two general situations can arise. First, specific features or details of the temporal or spatial behavior of the direct observables above may be of interest. Common cases of this type will consist of induction times or threshold behavior in kinetics or reactive flows. A second class of features could arise as functions or functionals of the raw system observables. An example in this category would consist of thermodynamic properties of a reacting flow. Given that such features are especially important in a laboratory

context, their sensitivities may be of more direct concern than that of the primary observables. With these sensitivities as a goal, there are several means of obtaining them. The choice of method, to some degree, also depends on how the system sensitivity equations are to be solved. In general, knowledge of the sensitivities 60i(r, t)/&xj(r’, t’) will always allow for subsequent calculation of the sensitivity of any characteristic feature in a system. However, if a small number of desired objectives are known beforehand, then an adjoint objective analysis may be preferred. Procedures i) and ii) fall in the former category and procedure iii) is the latter case.

i) Objective function jitting The objective function fitting approach is perhaps the most pedestrian of all the feature sensitivity analysis techniques, but it also has an important potential peripheral application to global sensitivity analysis (cf., section 3.2.2) which is not available by the other techniques. The basic idea behind the fitting approach is to first identify the features of interest which are assumed here to be characteristic parameters in the laboratory observable profiles as a function of temporal or spatial behavior. These parameters are then introduced into an explicit functional form and the parameters are adjusted to fit the calculated profiles. The sensitivity of the feature parameters may then be readily obtained. As an explicit illustration of this procedure, consider O(r, t, a) as the observation of interest where the parameter dependence is explicitly indicated. By an examination of the r and t dependence of this observation, it is assumed that meaningful characteristic features may be identified and an explicit functional form d(r, t, 8) chosen which contains the feature parameters &, &, . . . . By implication, the two forms of the observation are equivalent O(r, t,a)

= a(r, t,B(a)).

(12)

Eq. (12) implies a relationship between B and a.

H. Rabitz/ Local and global parametric ana&is of reacting flows

In practice, we will only know a solution of the model equations and the sensitivities at a reference point a0 in parameter space. This information will not be sufficient to determine the functional relation /3 = &a), however we may determine /3’ = &a’) at the system reference point and the corresponding sensitivity coefficients (@3Jaol,),0. In order to achieve this goal, the feature parameters in eq. (12) must be adjusted consistent with that relation. A technique of this type could employ minimization of the least squares functional

Minimization of R with respect to the feature parameters will yield the equation ~=JJdrdt[~(r,r,a)-B(r,t,8)] I

73

profile. Some special cases would be the location of concentration extrema, a threshold time, the width of a concentration excursion, etc. Such pointwise features may typically be associated with an explicit determining relationship 9 [W,

1, a),

sl = 0.

(14

Simple differentiation of eq. (14) would yield the desired feature sensitivity coefficients. As a concrete example, consider the location of a concentration extrema as a feature of interest in a temporal system. For a given concentration the extrema at time t * would be defined by the determining relation $(t,a)J

=O. t=t*

In order for this relation to be true, we must have Therefore, differentiation of the determining relation with respect to the system of parameters will produce the desired feature sensitivity coefficients at */hi = - (a*o/at aa,)/ ( a*O/at*) where it is understood that the derivatives on the right hand side are evaluated at the parameter reference point and the extrema location t*. It is evident that the feature sensitivity coefficient in this case may be determined from the original model equations as well as the sensitivity coefficients for the system. It is also apparent that second order sensitivity coefficients may also naturally arise when considering pointwise feature sensitivity analysis. These higher order derivatives can be calculated by procedures analogous to those discussed in 2.1. An extensive application of this technique has already been made in studying the moist combustion of carbon monoxide [8]. t * = t*(a).

x$$r,r,/+o.

(13)

The derivative in the integrand of eq. (13) may be explicitly evaluated by recalling that 6 has a known functional form with respect to its variables. Eq. (13) implies the existence of the relationship /3 = &a), but again, it must be recalled that O(r, t, a’) is assumed known only at the parameter reference point. Therefore, differentiation of eq. (13) with respect to one of the input parameters will yield an equation from which we may solve for the desired feature sensitivity coefficients ( J&/c?~~),o, In carrying out this last differentiation, it is evident that the system sensitivity coefficients aO,(r, t, a)/i3aj (or, if appropriate, their functional analog) will enter. The implementation of this overall procedure of feature sensitivity analysis is quite straightforward and in practice is only limited by one’s ingenuity in choosing simple but flexible functional forms d( r, t, #I). ii) Pointwise feature sensitivity analysis In many problems the features of interest can be indentified as characteristic points on a laboratory

iii) Adjoint feature sensitivity analysis

The pointwise special case of a would first consist as a functional of tions. The adjoint

technique above is, in fact, a more general procedure which of writing the identified feature the solution to the model equatechnique goes further and aims

14

H. Rabitz/

Local and global parametric analysis of reacting flows

to calculate the subsequent feature sensitivities without first calculating the standard system sensitivity coefficients. As an illustration, consider the following functional which is assumed to be a system feature of interest. 4t= I/F

[ 0( rrr, t”, a), r”, r”] dr”dr”.

(15)

For generality, we will assume that the system parameters are also distributed such that a functional variation of the objective is the desired quantity, g(

r, t) = //

“;:;;I

;;’

$

( rjr, t”) dr”dt”.

(16) The sensitivity entering into the integrand of eq. (16) may be identified by use of the explicit expression in eq. (11) to produce the relationship

where the reduced Green’s function g is given by

gw 0 = j-1dF( r”Jo’f “) x G( r”,

t”; r’, t’) dr”dt”.

(18)

(Note that the row vector g in eq. (17) is written as a column vector in eq. (18).) If an equation could be found for g without first obtaining the full Green’s function G, then a clear savings would be achieved by the calculation of a vector over a matrix and a reduction of the four arguments down to two. In order to achieve this savings note that the right hand side of eq. (18) is a projection upon the full Green’s function by a vector representing the particular objective of interest. Eq. (7) defining the Green’s function cannot be immediately utilized for this purpose due to the ordering of the arguments and the operators involved. However, the equation for the adjoint

Green’s function can be utilized. satisfies the differential equation G+(r’, t’; r, t)J+(r’, t’) = 16(r-

The adjoint

r’) 6(t - t’),

(19) where Gt( r’, t’; r, t) = G( r, t;

r’,

t’).

cw

The differential eq. (19) is in terms of the variables r’, t’ and therefore, the differential portion of J+ is understood to act to the left upon the adjoint Green’s function. Application of eq. (20) in eq. (18) along with the incorporation of eq. (19) will finally yield the desired differential equation for the reduced adjoint Green’s function g( r’, t’)J+( r’, t’) =

cw( r’, t’) ao

.

(21)

The computational advantages indicated above depend on two factors: (1) The objective must be known in advance of solving the problem; and (2) the number of objectives must be less than the number of direct observables arising in the original modelling equations. Furthermore, the advantage of this procedure becomes less attractive considering the ease of calculating sensitivity coefficients along the lines addressed in eqs. (4) and (5). Nevertheless, an adjoint objective function sensitivity analysis represents a powerful procedure of general utility.

3. Global sensitivity analysis The discussion in section 2 was confined to the consideration of local sensitivity issues in the sense that a system was being explored in the local vicinity of an operating point in parameter space. Such a local analysis is particularly useful for addressing the question, “What variables are important in a given system?’ Another class of problems is associated with the desire to know the

H, Rabitz/

Local and globalparametric analysis of reacting flows

solutions to a family of problems which are defined by the parameters spanning a wide region of the parameter space. Issues of this type naturally arise in the consideration of system optimization over control variables. In addition, for problems in which the parameters or some subset of them are entirely unknown, then exploration of a broad region of parameter space can be valuable. Furthermore, a simple understanding of the behavior at a given operating point in parameter space may be best achieved by exploring a broad neighborhood of that point. A further complicating issue arises from the fact that the nominally real system parameters should actually be treated as complex variables thereby raising the dimensionality of the space under consideration for exploration. Another class of problems of a global nature is related to the topic of objective function analysis. In this case, a region of parameter space is to be explored such that a constraint on the system is maintained. For example, in a hydrocarbon combustion system the fuel might be characterized by its carbon/ hydrogen ratio and it may be desired to know how values of this ratio influence various flame properties. Numerous other constrained searches could be envisioned. In general, global sensitivity analysis breaks into two broad categories: (1) matters of statistical error propagation; and (2) problems of global mapping. The somewhat different nature of these two problems will result in their separate discussion below. Before proceeding, it is important to realize that global problems of all types are inherently difficult and they will continue to be so categorized. This comment is especially true for global mapping. The reason for the comment is evident when it is realized that knowledge about a broad region of parameter space would accordingly be equivalent to knowledge about a broad class of problems. Such full knowledge represents a difficult but important objective. The scope of the problem can be further understood by realizing the traditional Monte Carlo or experimental design sampling of many points in parameter space would be exceedingly time consuming in spaces of high dimension. Therefore, the

15

general objective of global sensitivity techniques is to reduce the effort of these tasks to a practical level. Naturally, the complexity of global sensitivity issues has lead to the topic being less welI developed than the local techniques. In particular, the material below will focus exclusively on initial value problems which are in keeping with knowledge available at this time.

3.1. Global statistical metho& In general, global statistical methods aim to understand how statistical errors in the input of a model propagate through the dynamics and show up in the observable output or features. These statistical input errors could arise in the parameters of the differential equations as well as perhaps in the initial and/or boundary conditions when appropriate. Local sensitivity analysis can be applied to this issue, and it is a standard practice to use the sensitivity gradients for this purpose. Here we shall discuss more general techniques. 3.1.1. Statistical sensitivity analysis This procedure represents the most general statement of the statistical problem. A well posed problem would be specified by knowledge of the initial probability distribution function P(0, 0, a) of the observables and parameters and the goal is to calculate the probability distribution function P(t, 0, a) as a function of time [9]. This type of question is a standard one in statistical analysis. Given the problem defined by a system of ordinary differential equations

g=m4,

(22)

one may derive a diffusionless Fokker-Planck type equation for the desired probability distribution function

g+v’(fP)=o,

(23)

where the gradient operator has components vi =

H. Rnbitz / Local and globul parametric analysis of reuctingflows

76

a/m,.

Full knowledge about any statistical questions could be obtained by determining the distribution function from eq. (23). In particular, the average of any function of the observables or parameters g( 0, a) would be given by

(g(4)

=/I

dO,daP(t,

O,a)g(O,a).

(24

rewriting eq. (22) and expanding as indicted above g

=j((O)

+ At?, (a) + Aa)

Jf =f((Q,(a))+aAfl+zAa 1 + T &ABA0 [ + dAt2AaAaAa]

In many problems the observable-observable covariance matrix Co* or the observable-parameter covariance matrix Co” would be sufficient information.

coo=

((O-

C”“=((O-

(W(O-

(W),

(25)

@))(a-

(a))).

(26)

The partial differential equations in eq. (23) could be approached for solution by a number of means including direct explicit integration as well as the method of characteristics. Thus far, little work has proceeded along these lines, particularly for combustion-oriented problems.

3.1.2. Expected value analysis The expected value analysis (EVA) technique is closely related to that of statistical sensitivity analysis except that EVA attempts to calculate the mean observables and their covariance matrices without first obtaining the full probability distribution function [lo]. In general, this task cannot be achieved without making some approximation. The EVA approach consists of expanding the physical system differential equations about their nominal parameter vector as well as their nominal observation vector. In particular, these latter vectors are written as a = (a) + Aa and 0 = (0) + AO. In most problems the mean parameter vector (a) would be a known time independent quantity and the mean observation vector (0) would be sought after along with the system covariance matrices Coo, CO”. The analysis proceeds by

af + &AOAa + ... .

(27)

In practice, the expansion of eq. (27) must be truncated at some finite practical level with second order indicated above. Taking the expectation value of eq. (27) will eliminate the linear terms on the right-hand side and produce a differential equation for the trajectory of the mean observables coupled to the trajectories of the covariance matrices. Therefore, differential equations for these latter quantities must also be obtained and these may be established straight-forwardly by using the expansion in eq. (27) along with the definitions of the covariances in eqs. (25) and (26). Again, a hierarchy of coupled equations will be produced and the hierarchy must be closed in order for practical results to be achieved. In keeping with the present discussion, truncation through the covariance level will produce a closed set of coupled ordinary differential equations for the mean observables and their convariances along with that of the parameters. The physically important feature of this approach is exactly the coupling amongst the equations. In a local sensitivity analysis no such direct coupling occurs in the sensitivity equations and it is only introduced after their solution. In the present case, the coupling is an inherent part of the solution process. The price paid for this more complete information is a larger set of differential equations. In many cases, the increased size would not be a critical issue and calculations should be practical. An important issue to further explore is the accuracy of the results in relation to the hierarchical truncation procedure. Thusfar only limited calculations have been performed with this technique.

H. Rabitz/L.ocal and global parametric analysis of reactingflows

3.1.3. Fourier amplitude sensitivity test method This procedure has the same goals as that in parts 3.1.1 and 3.1.2, but the technique is quite different [ll]. The essential idea consists of first replacing each nominally constant parameter by a function of a new running variable s.

II

proach would simply consist of repeated calculations. It is naturally expected that any realistic modelling effort would entail some degree of this approach. However it is not expected to give a deep understanding of the mapping process. Beyond this technique, there are potentially a number of other procedures as indicated below.

a, + ajyj [sin (w/s)] . where wi, w2,. . . , would ideally be a set of incommensurate frequencies. The functions yj should be chosen as space filling curves such that as s varies, the parameters sample all regions of the parameter space. One may derive a differential equation for the functions yj given knowledge of the probability distribution function of the nominal system parameters. The result of the parameter replacement in eq. (28) is to introduce a new time-like variable into the output observables Oi(t) --, Oj( t, s). Statistically relevant sensitivity information is then achieved by Fourier transforming these new observables with respect to the variable s. The frequency spectrum is a direct reflection of the parameter dependencies including their nonlinear interactions. These latter effects would show up as overtones or combination frequencies associated with elements of the new parameter vector in eq. (28). In order to achieve a practical solution, the frequencies are typically chosen as being integral multiples of each other but spanning a wide frequency range. In this way the Fourier variable s need only be considered over a finite range, and the potential confusion between coincident overtones, fundamentals or combinations is still minimixed. Realistic calculations can become quite time consuming, since the original system of equations must be integrated numerous times corresponding to the number of discretized points of the variable s needed for the Fourier transform process.

3.2. Global parameter mapping As with the statistical methods discussed above, there are several possible approaches to global mapping. The obvious and most pedestrian ap-

3.2.1. Series continuation techniques Within the framework of local sensitivity analysis, the most direct global technique would attempt to take advantage of the availability of the sensitivity gradients. In particular, a Taylor series expansion could be considered.

O( a

+

ao

Aa) = O(u) + &da

i a20

+ZaaAaAa+

a**.

For sufficiently small variations this would presumably be an adequate approach, however for large scale variations there is no guarantee that the series actually converges. The use of eq. (29) clearly points out that a local analysis is best thought of in terms of its qualitative information about the question of the role of the system parameters. In contrast, a global analysis is concerned with the matter of the response to a finite change of the parameters. A possible approach in conjunction with eq. (29) for global mapping would involve re-expressing the series as a rational polynomial in the incremental parameter variation, Oj(a + Aa) = P;(Au)/Q~(Aa), where the numerator and denominator, respectively, are Nth and Mth order polynomials. Such Pad6 approximants have been successively applied in a wide variety of problems outside of the present context and they are known to provide a powerful means of continuation of series such as eq. (29) beyond their radius of convergence [12]. Some limited applications in kinetics have also shown promising results. The main difficulty with

78

H. Rabitz/ Local and global parametric

implementing this technique is the lack of availability of the necessary high order local sensitivity gradients needed to determine the coefficients in the Pad6 polynomials.

analysis of reactingflows

taken here, but for the sake of simplicity, we consider the desired objective to merely be a function of the original output observables. F= F[O(a,

3.2.2. Feature parameter approach to global scaling The technique of feature sensitivity analysis embodied by the relation in eq. (12) has an immediate spin-off application to global parameter scaling or mapping. Eq. (12), for our present purposes may be recast into the following form: O(r, t, a + Aa) = &(r, t, /3(a + da)).

t)].

(324

The simplest illustration of this technique would involve two parameters ai and a2 for variation with the remaining parameters fixed. The aim of this variation would be to find the relation a1 = ai that maintains the objective in eq. (32a). The maintainence of this objective implies that

(30) dF=~$~da,+$da,=O,

Again, as stated previously, this equivalence cannot be directly applied since we do not have full knowledge about the relation between /3 and a. However, the feature sensitivity analysis based on eq. (12) lead to knowledge of /3 and the sensitivities of /3 about the nominal operating point. Therefore, we may consider the expansion. p(a

+Aa) =/3(a) +XAa.

(31)

Substitution of eq. (31) into eq. (30) will yield a nonlinear scaling expression with respect to the parameters a in contrast to that achieved by eq. (29) using exactly the same first order information. This feature parameter scaling approach would be both computationally practical as well as likely to give acceptable results over an extended neighborhood of the system operating point. A clear example of this situation would arise in the singular perturbation problem of parameter dependence in oscillating flames. In those cases where the parameters influence the system frequency, a local sensitivity analysis will produce secular growth. In contrast, a feature analysis on the system frequency should be stable. 3.2.3. Objective mapping techniques The consideration of local objective sensitivity analysis in section 2.3 can basis of a global parameter space technique [13]. A functional approach

function form the mapping could be

I

’

from which equation

1

we may

2

obtain

(32b) the determining

(33) In the context of local sensitivity analysis, relations such as that in eq. (33) have also been referred to as derived sensitivity coefficients [2, 31. The relations may always be derived upon interchanging appropriate independent and dependent variables. In the present case a1 has become a dependent variable at the expense of constraining the observables through eq. (32a). In this framework global mapping is achieved by treating eq. (33) as a differential equation to be solved simultaneously with those determining the sensitivity coefficients with respect to ai and 0~~along with the original system equations. The result will be a trajectory in the space of ai and a2 along which the constraint in eq. (32a) will be maintained. As a final comment on this procedure, note that the expression in eq. (33) is time dependent through the various terms on the right-hand side. The derivative on the left-hand side was obtained by first calculating the observable sensitivities with respect to ai and a2 and effectively interchanging variables as discussed above. This procedure may be circumvented and direct differential equations may be obtained for &,/&x2.

H. Rabitz/

Local and globalparametric

3.24. A Lie approach to global mapping It appears that the most general approach to global mapping can be achieved by the consideration of Lie group techniques [14]. Lie group theory has long been used as a means for finding the independent solutions to differential equations. In the context of this paper, this usage would correspond to employing Lie theory for effectively solving the model equations. Here, we already assume that such solutions have been obtained by one means or another. Our goal here is to map the solutions throughout a finite region of the parameter space. This effectively corresponds to finding the solutions to many independent models. By confining ourselves again to initial value problems of the type described in eq. (22), the Lie approach may be understood in the following terms. We seek the transformations that perform (a, 0, t) + (a’, 0, t’)

(34)

such that in terms of the new variables we have the equation (35) The transformation eq.(34) has preserved the form of the original eq. (22). When considering the transformation in eq. (34) it is natural to think of the mapping a + a’ as causing a concommitant mapping of the solution 0 + O’, however the transformation t + t’ also, in general, needs to be considered. Such a dilation or contraction of time can naturally occur upon parameter changes. This would particularly be the case if a constrained mapping was involved. Therefore, the overall philosophy consists of seeking transformation operators that achieve eq. (34) while maintaining the validity of eq. (22) and (35). The task defined above is accomplished by tirst defining an infinitesimal differential operator

analysis of reactingflows

19

acting in the full space of all the variables. This operator is understood to produce the infinitesimal transformations corresponding to eq. (34) a’ = (1 +

daU)a,

0’ = (1 + daU) O’,

(37)

t’ = (1 + daU)t,

where a is referred to as the group parameter. Finite transformations would be generated by integrating the differential eqs. (37) or equivalently recognizing that the operator T = exp (au) would achieve the same task. At this point the coefficient functions entering into eq. (36) are unknown and must be determined. The first set of these functions 4 are indicated as special and only depend on the system parameter vector. This restriction is consistent with our desire to treat the transformation in eq. (34) as actually being brought about through the fundamental parameter transformation a + a’. Therefore, we restrict the vector I/J and furthermore, the explicit form of the elements in this vector is entirely at our disposal. These components dictate the nature of the transformation throughout the parameter subspace of the full problem. For example, if the choice JI=(l,l,l,...,) is made, then all the parameter space excursions are simple translations. Other choices could clearly be made if desired. At this point, it is necessary to determine the remaining coefficients l and +r. These are determined by demanding that the transformations leave eq. (22) invariant. In order to implement this demand, it is necessary to employ an extended generator since eq. (22) not only depends on a, 0, t but also on dO/dt. The extended generator can be shown to have the form

where

80

H. Rabitz/Local

and globalparumetrrc

Now demanding that

fi[g-_I(c9,a)] =o, will yield the sought after determining equation

(38)

The independent solutions to this set of partial differential equations will produce the corresponding independent infinitesimal generators for the transformations that preserve the original equations of motion. Several comments are in order concerning the structure of the determining equations in eq. (38). First, these equations are linear with respect to the unknowns [, 7. Secondly, there is a distinct relationship with the equations of sensitivity analysis. This connection may be understood by noting that two terms in eq. (38) may be combined into the total time derivative of c,,

Utilizing this identity in eq. (38) along with ignoring the time transformations r will lead to the following result:

The connection with sensitivity analysis may be understood by operating on eq. (22) with V= Z,$,a/&x, and making the connection Sj + Voj. When comparing these two results it must be recognized that the partial derivative residing in eq. (22) arises in recognition of the system pararneters being variables along with time t. Therefore,

ana!vsis of reactingflows

in essence the difference between the Lie approach and that of sensitivity analysis largely resides in how the operating equations are interpreted and utilized. In the case of traditional sensitivity analysis the sensitivity coefficients are sought in reference to a particular time dependent trajectory. On the other hand, the global Lie approach aims to map out a region of solution space. These differences of perspective are essential and translate into distinctions of methodology when treating the corresponding determining or sensitivity equations. There is thus far very little experience in the practical development of the Lie approach. It is also interesting to note that eq. (38) has a mathematical structure similar to that of eq. (23) and some of the same techniques may be useful for both problems. At this point only simple linear examples have been studied with the Lie parameter space mapping technique. However, even in these cases it is significant to note that examination of the infinitesimal generators gave valuable insight into the finite transformations without actually performing these transformations. In particular, numerical examples have shown how to continuously map from oscillatory to explosive regimes, and this latter transition was explicitly evident in the infinitesimal transformation generators. Two practical points about the Lie approach need consideration. The main point is that useful information can still be obtained without actually finding all the independent solutions to the determining equations. Lack of knowledge of all the transformations would imply that only certain types of mappings could not be achieved. Such a lack of full knowledge would probably only open up a subspace of the full problem for mapping. A second main point concerns a conjecture about obtaining approximate solutions to the determining equations. In practice, such approximations will be necessary since the equations will not likely be analytically soluble. Provided that any essential symmetries in a system are preserved, it is speculated that even relatively crude approximations to the determining equations may give adequate re-

H. Rabitz/

Local and global parametric analysis of reactingjlows

sults upon consideration of the system transformations. This conjecture is largely based on experience with Lie transformations outside the realm of the present context. Much basic developmental work needs to be done to make the Lie approach into a practical technique. Its primary attractiveness comes about through the observation that each question or objective of concern will generally require a new Lie generator. In other words, specific tools would be determined for each task at hand. This mode of operation should be contrasted with merely spreading sample points throughout the parameter space and resolving the problem. The latter technique, of course, works for all problems, but is optimal for none.

81

Acknowledgements

The author would like to express his appreciation to Geoffrey Ludford for inviting me to participate in the Special Year on Reacting Flows. His hospitality along with that of the chemistry department made the visit enjoyable and scientifically stimulating. Finally the author would like to acknowledge the Department of Energy, the Office of Naval Research, and the Air Force Office of Scientific Research for their support of various aspects of this research during the past several years.

References

4. Conclusion This paper, based on a series of lectures at Cornell University, aimed to summarize some particular and general aspects of sensitivity analysis of current interest in the kinetics and combustion communities. The actual lectures were accompanied by numerous numerical illustrations, many of which may be found in the available literature. Combustion problems which are inherently a complex interplay of chemical reactions and fluid mechanics provide an important area of application for sensitivity techniques. It is the author’s belief that much information concerning the important variables in reactive flows could be obtained by applying local sensitvity analysis methodology. This latter methodology, for some classes of problems, is already available for implementation. Global sensitivity issues are expected to receive increasing attention in terms of basic development and, ultimately for applications purposes. When considering all these matters it is significant to understand that the methodology is quite general and could be applied outside of the realms of combustion phenomena.

[l] R. Tomovic, and N. Vukobratovic, General Sensitivity Theory (American Elsevier, New York, 1972). [2] H. Rabitz, M. Kramer and D. Dacol, Ann. Rev. Phys. Chem. 34 (1983) 419. [3] H. Rabitz, in The Mathematics of Combustion, J. Buckmaster, ed. (SIAM, Philadelphia, 1985). [4] M. Smooke, J. Miller and R. Kee, Comb. Sci. and Techn. 34 (1983) 79. [5] W. Stewart and J. Sorensen, in Chemical Reactor Engineering, Proc. 6th European Symp., I-12,1976. H. Saito and L. Striven, J. Comp. Phys. 42 (1981) 53. B. Lojek, I.E.E.E. Proc. 129G (1982) 85. A. Dunker, Atomospheric Environment 15 (1981) 1155. Y. Reuven, M. Smooke and H. Rabitz, J. Comp. Phys., 1986, in press. [6] R. Hedges and H. Rabitz, J. Chem. Phys. 82 (1985) 3674. [7] M. Kramer, J. Calo and H. Rabitz, Appl. Math. Modelhng 5 (1981) 432. [8] R. Yetter, F. Dryer and H. Rabitz, Comb. and Flame 59 (1985) 107. [9] V. Costanza and J. Seinfeld, J. Chem. Phys. 74 (1981) 3852. [lo] Zi-ping Luo and H. Rabitz, to be published. [ll] R. Cukier, C. Fortuin, K. Schuler, A. Petschek and J. Schaibly, J. Chem. Phys. 59 (1973) 3873. [12] The PadC Approximant in Theoretical Physics, G. Baker and J. Gammel, eds. (Academic Press, New York, 1970). [13] M. Kubicek and V. Hlavacek, Numerical Solution of Nonlinear Boundary Value Problems with Applications (Prentice-Hall, New Jersey, 1983). [14] L. Hubbard, C. Wulfman and H. Rabitz, J. Phys. Chem., 1986, in press.