I. Linear Distributed Parameter Filtering: Observation Processes and Boundary Conditions for Engineering Systems. II. Statistical Analysis of Air Pollutant Observations and Model Predictions

Citation

Bencala, Kenneth Edward (1979) I. Linear Distributed Parameter Filtering: Observation Processes and Boundary Conditions for Engineering Systems. II. Statistical Analysis of Air Pollutant Observations and Model Predictions. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/r2fc-1237. https://resolver.caltech.edu/CaltechTHESIS:02172026-231434100

Abstract

Part I

The linear distributed parameter filtering problem is considered for systems of engineering interest. In particular, the boundary conditions and observation processes considered represent those encountered in practice. The intent is to consolidate the theoretical results relevant to engineering applications and to provide a unified structure to facilitate implementation of these results through a complete and consistent set of formal derivations. When necessary, previous results are extended.

The types of observation processes considered are

- continuous time, discrete space

- discrete time, discrete space.

The discrete time case includes both instantaneous and limited time average observations. And the discrete space case includes both pointwise and integral observations. Based on the orthogonal projection lemma, the minimum variance estimate is derived.

To examine the effect of stochastic boundary conditions, second order systems with mixed and Dirichlet boundary conditions are considered.

Part II

Aspects of air quality analysis that must be characterized with statistical methods are considered. The evaluation of air quality models is considered in detail . First, an evaluation framework is developed. The significance of evaluations of validity, accuracy and efficiency are discussed. Then, the specific aspect of accuracy assessment is addressed. The emphasis is placed on practical and objective methods. An extensive package of specific methods is set forth for use with air quality models. The package has been coded in FORTRAN. A description of the code is included.

In an attempt to increase the understanding of characterizations of long-term data, observed frequency distributions of air pollutant concentration levels are critically analyzed with respect to their statistical description. It is demonstrated that several common distributions can be used to fit observed data, one of which is the popular log-normal distribution. The observation that concentration distributions for all averaging times are approximately log-normal can be explained if the short averaging time data are themselves assumed to be log-normally distributed. The near log-normality of pollutant concentration frequency distributions can be explained on the basis of the near log-normality of wind speed distributions, although this explanation does not establish that wind speed distributions are solely responsible for observed concentration distributions. It is concluded that pollutant concentration frequency distributions are the result of complex phenomena and cannot be predicted exactly, but that the approximate log-normal character of the distributions is useful from a practical point of view and can be understood qualitatively on the basis of the relation between wind speed and concentration.

Item Type:

Thesis (Dissertation (Ph.D.))

Subject Keywords:

(Chemical Engineering)

Degree Grantor:

California Institute of Technology

Division:

Chemistry and Chemical Engineering

Major Option:

Chemical Engineering

Thesis Availability:

Public (worldwide access)

Research Advisor(s):

Seinfeld, John H.

Thesis Committee:

Unknown, Unknown

Defense Date:

23 October 1978

Funders:

Funding Agency	Grant Number
NSF	UNSPECIFIED
Environmental Protection Agency (EPA)	UNSPECIFIED

Record Number:

CaltechTHESIS:02172026-231434100

Persistent URL:

https://resolver.caltech.edu/CaltechTHESIS:02172026-231434100

DOI:

10.7907/r2fc-1237

Default Usage Policy:

No commercial reproduction, distribution, display or performance rights in this work are provided.

ID Code:

18372

Collection:

CaltechTHESIS

Deposited By:

Benjamin Perez

Deposited On:

23 Feb 2026 20:17

Last Modified:

23 Feb 2026 20:29

Thesis Files

PDF - Final Version
See Usage Policy.
34MB

Repository Staff Only: item control page

I. LINEAR DISTRIBUTED PARAMETER FILTERING : OBSERVATION PROCESSES AND BOUNDARY CONDITIONS FOR ENGINEERING SYSTEMS I I. STATISTICAL ANALYSIS OF AIR POLLUTANT OBSERVATIONS AND MODEL PREDICTIONS Thesis by Kenneth Edward Bencala In Partial Fulfillment of the Requirements for the Deqree of Doctor of Philosophy California Institute of Technology Pasadena, California 91125 1979 (Submitted October 23, 1978) ii. Acknowledgments I wish to thank my advisor, Professor John H. Seinfeld, for his advice and interest in my thesis work. The financial support for this research by the National Science Foundation and the Environmental Protection Agency is gratefully acknowledged. The list of persons at Caltech interested in systems theory and air pollution is far too long to write down here, and I would probably miss somebody. So, without particulars, I will merely thank all of those from whom I have learned while at Caltech. I also appreciate the assistance of the many members of the Caltech staff that I have dealt with during my stay here. In particular, the unselfish help and friendship of the following people has been very important to me: Edith Huang (Computing Center), Sally Ghan (Placement Office), Kathy Lewis {Chemical Engineering), Katie Foley, Rayma Harrison and Helen Lyons (Engineering Libraries) and, especially, a deeply felt acknowledgment to Lenore Kerner. For almost four years I shared an office with Tom Peterson and Sudarshan Kumar. I know not what pains and trials this arrangement may have meant for them, but for myself I feel I gained much. I thank them both. Many other students also contributed to the positive , sensitive and occasionally riotous spi rit among those of us in the Chemical and Environmental Engineering programs. Ajit Yoganathan, Gary Whatley, Vanis Yortsos, Freddy Gelbard and Dave Strand have each raised the qual ity of day-to-day life at Caltech. 1., 1., 1., . The incredibly undaunted faith of IT\Y parents, brothers, sisters and in-laws that I would someday graduate has been greatly appreciated. Before I finish this page , I must play the proud new father and acknowledge the sm.ile of IT\Y baby girl, Kristen. She is truly an immense joy. Finally, IT\Y best friend (and wife), Marcia, has been a constant source of support and counsel . I thank her dearly and most of all . iv. Abstract Part I The linear distributed parameter filtering problem is considered for systems of engineering interest. In particular, the boundary con- ditions and observation processes considered represent those encountered in practice. The intent is to consolidate the theoretical results rele- vant to engineering applications and to provide a unified structure to facilitate implementation of these results through a complete and consistent set of formal derivations. When necessary, previous results are extended. The types of observation processes considered are - continuous time, discrete space - discrete time, discrete space. The discrete time case includes both instantaneous and limited time average observations. And the discrete space case includes both pointwise and integral observations. Based on the orthogonal projection lemma , the minimum variance estimate is derived. To examine the effect of stochastic order systems with mixed and Dirichlet boundary conditions, second boundary conditions are considered. Part II Aspects of air quality analysis that must be characterized with statistical methods are considered. is considered in detail . The evaluation of air quality models First, an evaluation framework is developed. v. The significance of evaluations of validity, accuracy and efficiency are discussed. Then, the specific aspect of accuracy assessment is addressed . The emphasis is placed on practical and objective methods. An extensive package of specific methods is set forth for use with air quality models. The package has been coded in FORTRAN. A description of the code is included. In an attempt to increase the understanding of characterizations of long-term data, observed frequency distributions of air pollutant concentration levels are critically analyzed with respect to their statistical description. It is demonstrated that several common distributions can be used to fit observed data, one of which is the popular log-normal distribution. The observation that concentration distributions for all averaging times are approximately log-normal can be explained if the short averaging time data are themselves assumed to be log-normally distributed. The near log-normality of pollutant concentration frequency distributions can be explained on the basis of the near log-normality of wind speed distributions, although this explanation does not establish that wind speed distributions are solely responsible for observed concentration distributions. It is concluded that pollutant concentration frequency distributions are the result of complex phenomena and cannot be predicted exactly, but that the approximate log-normal character of the distributions is useful from a practical point of view and can be understood qualitatively on the basis of the relation between wind speed and concentration. vi . Table of Contents ACKNOWLEDGMENTS ii ABSTRACT . . . . iv PART I. LINEAR DISTRIBUTED PARAMETER FILTERING: OBSERVATION PROCESSES AND BOUNDARY CONDITIONS FOR ENGINEERING SYSTEMS . . . . . CHAPTER 1: 1 INTRODUCTION. 2 ...................... 1.1 Introduction 1. 2 Distributed Parameter Filters : Summary. 1.3 Specification of the System and Observation Processes . . . . . . . . . . . . . . . . . . 20 CHAPTER 2: 5 DISTRIBUTED PARAMETER FILTERS: DETERMINISTIC BOUNDARY CONDITIONS. 2.1 Introduction . . . . . . . . 2.2 Specification of the Filters: and Criterion . . . . . . 2.3 Summary of Results 2.4 Derivation of Filters • • 29 30 Notation • • 31 . . . . . . . . . . . . . . . . . . . . . . . . 33 37 ............ 2. 4.1 Continuous Time Observations 2. 4.. 2 Discrete-Instantaneous Time Observations 2.4.3 Discrete-Average Time Observations CHAPTER 3: .3 37 • • 44 ......... • • 49 DISTRIBUTED PARAMETER FILTERS: STOCHASTIC BOUNDARY CONDITIONS 3. 1 Introduction . . . 3.2 Summary of Resul ts . .... . . . . . . . . . . • • • • • 5.3 • • 54 . . . . . . . . . . . . . . . . . . . 55 vii. Table of Contents (Continued) 3.3 Derivation of Filters . . . . . . 3.3.1 Mixed Boundary Conditions . . . ... • • • • • 58 ...... 3.3.2 Dirichlet Boundary Conditions. • • • • 63 . . . . • • • • • 67 CHAPTER 4: SUMMARY . . . . . . . . . REFERENCES ...................... • • 76 • • • • 78 .... .. APPENDIX B: CERTAIN LUMPED PARAMETER FILTERS . . . . . . . . . . Introduction . . . . . . . . . . . . . . B.1 ... Discrete-Instantaneous Observations . . . . . . B.2 APPENDIX A: ORTHOGONAL PROJECTION LEMMA . . • • 82 • • 86 • • 87 • • 91 8.3 Discrete-Average Observations . . . . . • • 94 B.4 Continuous Observations with a Deterministic Input and Non-Zero Mean Initial Condition . . . . . . . . 95 PART II. STATISTICAL ANALYSIS OF AIR POLLUTANT OBSERVATIONS AND MODEL PREDICTIONS . 99 CHAPTER 1: INTRODUCTION .. 100 CHAPTER 2: EVALUATION OF AIR QUALITY MODELS: THE ACCURACY OF PREDICTIONS RELATIVE TO OBSERVATIONS. 102 Introduction 103 Evaluation Framework. 105 viii . Table of Contents (Continued) Accuracy Assessment. 10.9 Introduction . . . 109 Assessment Levels 1 and 2 . 11 P, Summary . . References APPENDIX. CHAPTER 3: AQMAAP, Air Quality Model Accuracy Assessment Package . . . . . . . . . . . . . 128 ON FREQUENCY DISTRIBUTIONS OF AIR POLLUTANT CONCENTRATIONS . 144 Introduction . . 145 Representation of Concentration Data . 146 Effect of Averaging Time on Frequency Di stri buti ons . . 147 Influence of Wind Speed on Instantaneous Concentrations . . . . . . . . . . 150 Influence of Wind Speed ·and Mixing Depth on Mean Concentrations . . 151 Summary . . • 152 Conclusions. • 15.3 References • 1 5.3 APPENDIX. The DistributiDn of the Average of Two Correlated Loq-normal Variates . . . . . . . . 15?, 1. PART I. LINEAR DISTRIBUTED PARAMETER FILTERING : OBSERVATION PROCESSES AND BOUNDARY CONDITIONS FOR ENGINEERING SYSTEMS Chapter 1 Introduction 3. Research concerned with the linear distributed parameter filtering problem has matured to the stage at which a unifying presentation of the problems, the derivations and the results is beneficial. As our orientation is to eventual applications, we emphasize completely specified results for the problems we consider. We define the linear distributed parameter filtering problem in a consistent manner for a variety of observation processes and boundary conditions. All of the filter derivations are based on the orthogonal projection lemma. In this chapter we introduce the filtering problem and specific problems and aspects to be studied. 1.1 Introduction In this work we develop linear distributed parameter filters for systems with boundary conditions and observation processes of the types encountered in the analysis of engineering systems. This first section consists of a general introduction to the filtering problem, with particular reference to alternative types of boundary conditions and observation processes. ' The purpose of such a general introduction is to provide a basis for the summary of previous theoretical work presented in Section 1.2. The particular systems a~d observation processes we consider are specified in the final section of the chapter. The linear distributed parameter filt~ring -problem deals with I ,, estimating the state of a stochastic syst~m ~overned by an equation of the type 4. where u is the state vector of the system, w a stochastic input, Lx a linear partial differential operator, and Ga matrix. be more carefully defined later.) on observations of the system. (The system will We desire to estimate the state based Due to the fundamentally random nature of physical processes, these observations are corrupted by random errors. The observations are expressed in the form Z =Hu+ v where z is the observation vector, v a random observation error, and H a bounded linear operator. There are a number of possible variations to this basic problem. We are concerned, in particular, with the different types of boundary conditions and observation processes. In the derivation of the dis- tributed parameter filters, it is frequently assumed that the boundary conditions are deterministic, i.e., 8 u = h X b where 8x is an appropriate linear operator and hb is a known deterministic function. A more qeneral form of the boundary conditions includes both detenninistic and stochastic inhomogeneities, 8xu =' hb + wb where wb is~ stochastic functiori~ Although there are many systems for which deterministic botJndq.ry conditions, apply, for the majority of I real sy'stems stochastic boundary· inputs w.ill be present . 5. Relative to the observation processes, the most frequently considered process in theoretical developments is, in fact, physically impossible to implement. In this case it is assumed that observations are taken continuously in time and at every point in the spatial domain In practice, observations are generally made at discrete locations, with the sampling being point-wise or integral. For certain state variables it is possible to make observations continuously in time, but frequently it is preferable, or necessary, to sample discretely in time on either an instantaneous or an average basis. 1.2 Distributed Parameter Filters: SuJTmary The last decade has been a period of active research on the distributed parameter filtering problem. Through both rigorous and formal approaches a reasonably complete theory and set of algorithms have been developed. In the following tables this progress is summarized. Our focus here is on the formal development of filters useful for linear engineering systems. For completeness, references to rigorous theories and to filters for nonlinear, systems are provided in Tables 1.6 and 1.7 respectively. Some background references to the study of distributed parameter systems are given in Table 1.1 , A surm1ary of the references to formal approaches for the linear filtering problem is given in Table 1.2. Tables 1.3 - 1.5 these references ~re further discussed. In 6. Table 1. 1 Distributed Parameter Systems Background References Reference Comments Athans 1970 Discussion of tentative approaches to unsolved problems in the practical implementation of the theories of distributed parameter systems. Robinson 1971 Survey of control--theory and applications. Curtain 1975a Review of filtering theory. rigorous theories. Polis and Goodson 1976 Survey of parameter identification--theory and applications. Ray 1977a Survey of filtering--t.heory and applications. Ray 1977b Survey of control, filtering and identification-applications. Tzafestas 1978 Extensive review of all previous· results in distributed parameter filtering. Focuses on the ?. Table 1.2 Fonnal Linear Distributed Parameter Filters Classified by the Type of Observation Process Considered OBSERVATIONS IN TIME Discrete-Instantaneous Continuous (Table 1.3) Tzafestas and Nightingale 1968a, b Tzafestas 1969, 1972a, b, 1973 Meditch 1971 Atre 1972 :ontinuous Atre and Lamba 1972a, b, c* Kumar and Sage 1972 Shukla and Srinath 1972 Padmanabhan and Colantuoni 1974 Reddy and Rajamani 1975a, b Lee 1976 OBSERVATIONS IN SPACE (Table 1.4) Thau 1969 Meditch 1970**, 1971 Sakawa 1972 Discrete Shukla and Srinath 1972 Tzafestas 1972a Padmanabhan and Colantuoni 1974 *Observations continuous on the boundary. **A single integral observation. Padmanabhan and Colantuoni 1974 (Table 1.5) Thau 1969 Padmanabhan and Colantuoni 1974 Aidarous and Ghonaimy 1976 Stochastic Deterministic Shukla and Srinath 1972 Stochastic Deterministic Kumar and Sage 1972 Atre and // Lamba 1972a Atre and Lamba 1972b Atre and Lamba 1972c Atre 1972 Meditch 1971 Reference Tzafestas and Nightingale 1968a Tzafestas and Nightingale 1968b Tzafestas 1969 Minimize the error covariance via the innovations approach. Minimize the error covariance via the innovations approach. Introduce the boundary Green's kernel to deal with the boundary conditions. (Observations continuous on the boundary.) Minimize the error covariance via the innovations approach. The boundary conditions are given as part of the problem specificaUon. Minimize the error covariance via the orthogonal projection lemma. Linear Distributed Parameter Filters-Continuous Time and Continuous Space Observations Boundary Conditions Aeeroach Deterministic Minimize the error covariance via the orthogonal projection lerrma. Deterministic Find the conditional expectation of the state via characteristic functionals. Stochastic Introduce the boundary kernel matrix to deal with the boundary conditions. Solve a deterministic optimal control problem via the De:t:erministic sweep method. Stochastic Introduce extended operators to deal with the boundary conditions. Deterministic Minimize the error covariance via the maximum principle. c:o Determ,i ni sti c Reddy and Rajamani 1975a Reddy and Rajamani _1975b Lee 197)5 / / Stochastic Padmanabhan and Colantuoni 1974 ' Deterministic Tzafestas 1973 Deterministic Determini$tic Boundary Conditions Deterministic Deterministic Reference Tzafestas 1972a Tzafestas 1972b Table 1.3 (continued) Minimize the error covariance via the Martingale representation. Ae_eroach Maximize the marginal density functional Extend previous results to include colored observation noise. Minimize the error covariance via an extended GaussMarkov Theorem and the innovations approach. Take the limit of the filter for discrete time observations. Introduce the Green's function to deal with the boundary conditions. Extend previous results to the case of hyperbolic systems with colored observation noise. Extend previous results to include time delays. '° Deterministic Deterministic Meditch 1970 Meditch 1971 Stochastic Deterministic Deterministic Stochastic Sakawa 1972 Shukla ,,and ·~ _.· Sri nath 1972 Tzafestas 1972a Padmanabhan and Colantuoni 1974 .. Boundary Conditions Deterministic Reference Thau 1969 Space Observations Ae_e.roach Minimize the error covariance via the calculus of variations. Solve a deterministic optimal control problem via the sweep method. (A single integral observation.) Start with the filter for continuous space observations and specify an assumed weighting of the observations. Minimize the error covariance via the calculus of variations. Introduce a generalized function to deal with the boundary conditions. Start with the filter for continuous space observations and specify an assumed form for the discrete observations. Start with the filter for continuous space observations and specify an assumed form for the observation noise covariance . Take the limit of the filter for discrete time observations. Introduce the Green's function to deal with the boundary conditions. Continuous Time and Discrete Linear Distributed Parameter Filters-- Table 1.4 '---' c::, '; Reference Thau 1969 Padmanabhan and Colantuoni 1974 Aidarous and Ghonaimy 1976 Deterministic Boundary Conditions Deterministic Stochastic A2..eroach Minimize the error covariance via direct differentiation. Find the conditional mean. Introduce the Green 1 s function to deal with the boundary conditions. Find the optimal gain in tenns of orthonormal basis functions. Discrete-Instantaneous Time and Discrete Space Observations Linear Distributed Parameter Filters-- Table 1.5 N N ! ,,,,,. Kushner, 1970 Falb, 1967 Balakrishnan and Lions, 1967 Reference •u(x', t)d0x x. t) dt +G(x, t)dw 1 du(x, t) = [Lxu + fo L(x bounded operator. LX t)dw 2nd order 1 inear spatial operator, + G(x, fr. du(x, t) = [Lxu - h(x, t)Jdt I. A: It an infinitesimal generator of a strongly continuous semigroup. du= A(t)udt + G(t)dw A: (t) + Au(t) = 0 * St,stem dab= Ba(t)abdt + Gab(t)dwa :v : cononnal derivative + Gb(x, t)ab(t) :~ + B(x, t)u = hb(x, t) u(x, t) and Lxu(x, t) + 0 as x + ao Boundart, Conditions Linear D1str1buted Parameter Filters - Rigorous Approaches Table 1.6 dy = Jolt + 1t Gv(t)dv + M(x, t')u'(x, t')dox}t' [fa0M(x, t)u(x, t)daox}t dv y(t) = dz= M(t)udt + Gv(t)dv z(t) = M[u(t)J + v(t) Observations h..l C:--:0 13alakrishnan, 1974a / / Phillipson, 1971 ~~ {t) + A(t)u(t) = h(t) + G(t)w(t) Bensoussan, 1971 ;f (x, t) = Lxu(x, t) + h(x, t) u(x, t) = ab(x, t) I. Continuous time/continuous space z(x, t) = u(x, t) + v(x, t) Observatiens . M: linear bounded transfonnation continuous semigroup. A: infi nitesimal generator of a strongly Also inexact observations of u(x, o) and ab are taken. z(t; w) = Mu(t; w) + Gvw(t; w) ~v: cononnal derivative No apriori information available about ab. = ab(x, t) au Continuous time/discrete space II. av (x, t) = ab(x, t) II. z(xi' t) = u(xi, t) + v(xi, t) au No state observation III. av (x, t) + B(x, t)u6<, t) III. z(x, t) = O I. Boundari'_ Conditions ~ (t, w) = Ax(t; w) + Gw(t; w) II. - a2u (x, t) = L u(x, t) + h(x, t) at 2.• X 1,c: 2nd order linear spatial operator. I. A: linear operator, allowing for a wide class of distributed parameter systems. y{O)=yo+wo s1.stem Reference Table 1.6 (continued) f--l ~ , • Kumar ani:I' Seinfeld, 1978 I Koda, 1975 z(t) = M(t)u(t) + v(t) M: Bound linear operator. Specified to yield spatially-independent integral observations and discrete point observations. du _ A ( ) ( •) ( ) ·ctt - u t + G t wt A: infinitesimal generator of a strongly continuous semigroup. M: Continuous mapping of a Hilbert Space into the measurement space. A: infinitesimal generator of a . ~-- semigroup. + Gb(t)dv(t, w} dz(t, w) = M(t)u(t, w)dt M: Bounded linear transformation m~pping L2(D) into another Hilbert Space. z(t) = Mxu(x, t) + v(t) Observations z(tk) = M(tk)u(tk) + v(tk) unbounded operator which generates an evolution optrator. u(x, t) = wb(t) Boundary Conditions ~ = Au(t) + G(t)w(t) A(t): + G(t)dw(t, w) du(t, w) = A(t)u(t, w)dt Curtain, 1975b at v2: Laplacian ~ = v2u System Balakrishnan, 1974b Reference Table 1.6 (continued) N ~ Lamont and Kumar, 1972 Hwang, Seinfeld and Gavalas, J972 and Hwang, 1971 Seinfeld. Gavahs Nightingale, 1969 Tzafestas and Reference System ~-4)• w(x , t) ax au / u ) N t, x, u, "ax•~•a + w(x, t) • ax ( X N : nonltnear second order spattal operator. MLl1. . NX (u, x, t) + w(x, t) at 0 Ha(t, a) + w (t) *• -ar---'" • au(x, t) // ~ - H(t, x, u, Wx: nonlinear spattal operat? ~ ----rr- • N,(u, x, t) + G(x, t)w(x, t) au(x. t) X *) • wb(t) xidO , I I nonl i near spatial operator. ,cl! z(x , t) • M(u, x, t) + v(x, t) BN : x X Continuous t1me/Cont1nuous space Stoctiast1c 01screte space also treated. xc[O, II z(x, t) • H(t, x , u) + v(x, t) Continuous t1me/Cont1nuous space f. BN (u, x, t) • wb(x, t) ift • Ba(t, ab) + w,b(t) dab e(t, 'b' u, General with stochastic inputs. cit• B0 (t, ab) + w,b(t) xcO uH(t) • Lu T(x, t), ... , uT(xp' t) ]T z(t) • H(t, uM) + v(t) dab Continuous ttme/Otscrete space z(x, t) • Mxu(x, t} + v(x. t) I. ab' u, 'ax au) • 0 B\t• ~~~! ~~:~ spat ta l Observa t 1ons Cont tnuous t tme/Cont tnuous space General w1th stochastic tnputs. BNX: O<!termlnlst1c RN (u, x, t) • 0 Boundary Conditions Nonlinear Distributed Parameter Filters Table 1.7 AEP_roach Solve a determ1n1sttc control problem vh 1nvar1ant 1mbedd1ng. by the ad hoc specification of the observation wt f ght 1ng ma tr1x. Discrete space observations treat ed Solve a determtn1sttc control problem v1a 1nvartant 1mbedd1ng. and take the 1 tm1 t back to the cont 1nuous case. apply the 1"11ped parameter filter D1scret1ze the systl!fl'\ spatially , Max tint ze the l f kel 1hood functional via dynamfc progra11Tnfng. Coornents The parameters a and ab are also estimated. Includes smoothing solut lon. 1:\~!~a~~~~:kd. Includes smooth1ng solution . N (Ji Yu , Setnfe ld and Ra y , 1974 i% (x, t) • N,(u, , , t) • w(,, t) ,/ . .,t 1 . -·- /Q •2 1 -t wa( t) l•• u(x , t), . .• , u(xn , t),t ) r N (u, x' • t)d x ,. l a• • Na IT • w( ,, t) i"f (x, t) • - N1(,, t) }J. .. N2(u, x, t ) N,., : noollnr:ar spa tial operator . Sher ry and Shen, 1973 second order spdt1al oper ator. ~ (x, t) • N,l u, x , t) • G(x, t)w(x, t) N1,,: ~ • N,(u, x , t) • G(x, l)w( x , t) Si_stem Tzafestas. 19 72 a Sahg a l and Webb• 1972 Reference Table l.7 (continued) X ftrst order spat 1al operator. X o~ rat or . non li near spatial u(O, t) • B lx( t )) .11 eN : Oe t enutnt s ttc eN (u, ,, t) • 0 X Oetenn tntstt c eN (u, x, t) • 0 x bN : 8N (u, x, t) '"0 Boundari_ Cond it t on Hanogene:ous Observa t1 ons 't•D xi:D xcO • v(t) x £(0 , 11 1 l • loMiu(x', t ), x ' , t)dx ' z(t) • H (a(t), u(, , t), ... ,u(xp' t)) 1 1 space Cont fnuous ttme /Otscre te and Integral • v(x t, tJ) z( , 1 , t j ) • M( x 1, tj)u( x 1 , tj ) Discrete t1me/0 1sc ret e space z( x , t) • M,(u, x , t) • v(x, t ) Con t 1nuou s t 1me/Cont t nuous space z(x, t) • M(u, x , t) • v(x, t) Cont1nuou s t1me/Cont1nu ous space A.ep_roach employing d1fferent1al sens 1ttv1tie s. S01ve a detenn1ntsttc control pr oD l em Solve a detennlntstlc control problem vt• the sweep method . estt.,.te and app ly the l inea r f11ter, L 1near1ze the sys tem around the ft 1ter So lv e a detenn1n1st1c con tro l problem v1a tnvartant 1mbeddtng. Clllments also estimated . The pdrdn~ l t! r a 1s solution. Includes smooth i ng solution. Includes smootnfng N C') Ray and Setnfeld, 1975 Seinfeld , 1975 Reference Ajtnkya, Ray, Yu, and Table 1.7 (continued) OU Sys tem Niu. ~ . ~ 1 0 0 x, t \ dx + w(x, t) du 2 u ' x, 'aa7°, t) J'"l,2, .... r+l.aj-l<x<~j + wJ( x, t) • Hj ( uj,, ? . wa (x, t) n1 (x, t) t t) H.(,, u(x , t). ... ,u(x 11 , t), t) 1 *• if (x, Boundary Cond t t tans 0 .:aQ wb, r t 1 "' O •u r + l --rx-• t) x • ar + l .. l t 8 r + l ( ur t l' xa a Bo ( ul" !x• •ul t) + wbO • 0 *•'• t) • wb(t) Stochastic s(u, General w1th stochasttc tnputs. Observat 1ons prising the moving boundaries and state variables. respectively. a and u denote the vectors com- z(t) • H(a, u(x , t). .• .,u(xp' t)) 1 • v( t) •1•10, 1} Cont tnuous t tme/Otscrete space Discrete ttme also treated . + v( t) z(t) • H(a(t). u(xtt), ... ,u(xp' t), t ) Continuous ttme/Otscrete space Approach Solve a detennfn1st1c control problem employing d1fferent1a1 sens1t1vtt1es. Discrete time observations treated by the ad hoc spec1f1cat1on of the observat1on we1ght1ng matrix. Solve a detenntntsttc control problem employing dtfferentlal sens1tlvlttes. Ccmnents to case of moYtng bound a r 1es. An eAtenston of the prect:d1ng two papers also estt11dted. The parameter a 1s ~ N r Tufestds, 1976 Ray c1.nc.J Sc1f1feld. 1975 (con 1 t} Rtftreoce S1.s ten1 J • l,. , ., r ,/ *-:.<(•, t)dt t G•(a, t)dw,(t) t (\(u, a, ~ • t )dw(x, t) du(x, t) • N,(u, •• , , t)dt t waj ( l) ".."J. a, (aJ' t), ~ a, (,J' t), •J' t ) ~ (t) • N.(u/•J• t). uj t 1(aj' t), Tht locat ions. aj, of tht r mo vtng bouudartes are determined by Tdble I. 7 (continued) au t t) wb J. • 0 ,--t,!---1- d•b t G,b(•o· t)dw.b( t) at• B•(•b • t)dt X BN (at,, x , t) • wb(t) General w1th stochastic tnputs. • J, l) t wbj • 0 j • I, . . . , au . au _J_ • -t l 8J' ( "J' "J t l' ~, aj , 8i ( uj, uj + l' i-?- , ~ au Bourldarz Cond tt1 ons G(,, •• t)dv(x, t) xcD Discrete space 1s also treated. t dz(x, t) • M(u, •• ab' z, x , t)dt Cont 1nuous t lme /Con t 1nuous sp ace Observa t ton s vation we1ght1ng matrix . tne ad hoc spec 1f1cat1on of the obser - Dtscrete space observati ons treated by equat ton. state via the fokker-Planck-Kolmogorov Ftnd the condtttonal expectton of the A.e.eroach ::~ i!Rt~~ 'i so The parameters a Coornents '-> a::, 1976 TZofostas and N1ghtlngale. Reference Srstem t G.11.(u, x, t )dw(x, t) ~ (x, t) • N,(u , x, t) d t Table 1.7 (continued) Uounda r d9nQ t~ ton s Oeterm1,dstfc t!x{u. x, t) • 0 Observdt I ons u ). proce ss are a fun ct 1on of the state, (the Stdt1st1cs of th~ count1ng Doubly stocna sttc counting process Continuous t1me/01screte space AE£roach Kolmogorov equation. the state v1a the Fok.ker-Planck- Find the cond1t1onal expectation of Coomen ts >---., '° 20 . 1. 3 Specification of the System and Observation Processes As stated in Section 1.1 we are interested in the linear distributed parameter filtering problem for systems with stochastic boundary conditions and for observation crete in space. processes that are dis- In this work we present a unified discussion of these particular problems. First, in Chapter 2, we treat the problem of discrete spatial observation processes for linear systems with deterministic boundary conditions . Three types of discrete spatial observation processes are considered . They are : - continuous time - discrete-instantaneous time - discrete-average time . Then, in Chapter 3, we show how these results can be extended in specific cases to systems with stochastic have chosen boundary conditions. We to consider a specific class of l i near systems and two types of bounda ry conditions. specify the results. This allows us to completely In most previous· work , consideration of more general cases has necessitated leaving unspecified functions in the fina 1 result . Tables 1. 8 and 1. 9 summar i ze the cas es considered in t his work and 1 list other references , if any, also contaj ni ng the f i lters . 4.1. Table 1.8 Linear Distributed Parameter Filters Considered in Chapter 2 system: ~~ = Lxu + G(x, t)w(x, t) + h(x, t) Boundary Conditions: Deterministic, sxu = hb X£3D References Also Containing Filter Observation Process Discrete Space Continuous Time 0 < t x£0 z(t) = JD M(x, t)u(x, t)dDx + v(t) As indicated in Section 1.2, this filter is widely available. (For example, see Sakawa, 1972 or Padmanabhan and Colantuoni, 1974.) Discrete-Instantaneous Time Padmanabhan and Colantuoni,1974. f z(tk) = DM(x, tk)u(x, tk)dDx + v(tk) Discrete-Average Time j z(tk) = ftk M(x, t)u(x, t)dDxdt tk - 1 D None Table 1.9 Linear Distributed Parameter Filters Considered in Chapter 3 2 au = a ( x, t) -a u + a (x, t) ax+ au a (x, t)u + G(x, t)w(x, t) + h(x, t) af 2 1 0 2 System: ax 0 < t 0 < X < 1 Boundary Condition References Also Containing Filter Mixed [ ¢1(t)a 2(1, t) ~~ + e1(t)u] X = 1 None for the specific system considered. Another second order system has been considered by Sakawa (1972). Dirichlet e1(t)u = h (t) + w (t) X =1 1 1 ; None e1 nonsingular Observation Process: Discrete Space; Discrete and Continuous Time We now specify the system and observation processes considered. Then-dimensional state vector, u(x, t), is governed by the partial differential equation G(x, t) w(x, t) + h(x, t) defined fort> 0, XED. (1.1) The domain, D ,is a connected subset of!- dimensional Euclidean space Rl with boundary surface ao. Lx is a linea~ partial differential operator with respect to the spatial coordinate x. w(x, t) is an n-dimensional stochastic input, white in time, his a known n-dimensional input, and G(x, t) is a known n x n matrix. The initial condition for (1.1) is (1.2) and the boundary condition * is (1.3) where Bx is a linear operator . When the stochastic terms are absent, the problem is assumed to be well-posed. *As indicated in Table 1.9, only for a specific class of these systems is the full stochastic boundary condition problem considered. 24 . The initial condition is unknown, with its first two moments given, ( 1. 4) The stochastic inputs have the following properties E{w(x, t)} = 0 (1. 5) E{w(x, t) wT(x', t')} = Q(x, x', t) o(t - t') (1. 6) where Q(x, x', t) and Qb(x, x', t) are non-negative definite n x n matrix functions. (The matrix S(x, x', t) is ' said to be non-negative definite, in a generalized sense, if ff nT(x")S(x", x"', t)n(x"' )dDxudDx"' for all n(x).) DD ?: 0 We are interested in observation pro~esses in which there are p discrete observation locations. z 1(t) cl'~~ot'e~ them-dimensional observation vector at the ; th observation location. Observations may be made in the domain or on the boundary. We employ the compact notation (Kumar 25. and Seinfeld, 1978), todenote the overall pm-dimensional observation vector. The three types ofobservations * considered are continuous time z(t) = f M(x', t) u(x', t) dDx' + v(t) (1. 7) 0 discrete-instantaneous time (1.8) discrete-average time** f M(x', t') u(x', t') dDx,dt' + v(tkl (1.9) D M(x, t) is the known overall pm x n observation matrix, and it represents a matrix consisting of the p observation m~trices at the individual locations, i.e., The following chart explains the temporal observation process nomen- clature. 26. Temporal Observation Process Nomenclature System Continuous u(t) Observation Sampling Process Transmission Continuous z(t) z(t) The recording of observations may be considered as a two-step process; sampling the system and transmitting the observations. For a continuous system, we may sample a system either in an instantaneous or in an average manner and then transmit the result either continuously or discretely. We consider three of these processes. The fourth, continuous-average or moving average as it is commonly referred to, is not considered. **in the case of discrete-average observations the observation noise includes both sampling and transmission noise . The noise occurring in the sampling process is averaged over the interval tk to tk -l · Thus, Equation (1.9) may be expressed tk tk vS(t )dt + vT(tk) z(tk) = J J M(x' ,t')u(x' ,t')dDx,dt' + f tk-1 D 1 tk-1 1 ( 1. 9a) where v(tk) = tk J vS(t')dt' + vT(tk)' tk-1 vs is the sampling noise and vT is the transmission noise. defined as independent white noises with properties E{v 5(t)v~(t')} = Rs(t)o(t - t') For vs and vT and E{vT(tk)vi(tl)} = RT(t)okl ( 1.12a) 27 . The observation noise, v, is a pm-vector, assumed to be zero-mean and white. For continuous time observations the properties of v(t) are E{v(t)} = 0 E{v(t)vT(t')} = R(t)o(t - t') = R1, 1(t)•·•R 1, p(t)] : : o(t - t') (1.10) [R 00 p, 1 (t) ■ R p, p(t) where Risa known positive definite pm x pm matrix and (1.11) Similarly, for discrete t1me observations the properties of v(tk) are (1.12,) where Risa known positive definit'e pm x pm matrix and i E{vi(tk)v}(t 1)} = Rij(tk ~okl" (1.13) Finally, u0 , w, wb and v are all defined ,to -be independent of each other . The observation processes have been expressed as integrals over the 28 . domain D. In applications one usually has point-wise observations. By taking the limit to small volumes of integration we can represent such a process. The observation matrices become i = 1, ... , p and the observations are i = 1, ... , p. 29 . Chapter 2 Distributed Parameter Filters: Deterministic Boundary Conditions 30 . In this chapter we present the linear distributed parameter filters for the three types of discrete space observation processes defined in Chapter 1. We specify the filtering problem, summarize the results, and derive the filters. 2.1 Introduction The system and observation processes we consider are as given in Section 1. 3 with detenn-1nistic boundary conditions. Stochastic boundary conditions are considered in Chapter 3. In Section 2.2 we introduce the necessary notation and define the estimate we seek. The filters for the three observation processes are presented in Section 2. 3 and derived in Section 2.4. The derivations follow the formal pro- cedure of Tzafestas and Nightingale (1968a) and are based on the orthogonal projection lenma (see Appendix A) . There are three differences in the exact specification of the problems considered in the derivations and in the summary, Table 2.1. These differences arise because of differences in the p4rposes of the der i vations and the summary. In the derivations we wish to present the essentials of the filtering probiem uncluttered by slight extensions. However , in the summary we wish to-present the filters in forms more useful for applications. following table. ' The three differences are stated in the i ' 31 . Derivations SuITT11ar~ Spatial Observations Integral Point-Wise Detenninistic Inputs h = 0 hb = 0 h hb uo = 0 uo Mean of Initial Condition 2.2 Specification of the Filters: Notation and Criterion We begin by defining the following tenns: u(x, tit') "'u(x, t) = u(x, tit') = u(x, P(x, tit'; x', t"lt"') = P(x, X', t It I) = P(x, P(x, XI, t) = w( x, t It 1 ) "' u(x, tit) "' t) - u(x, tit') Optimal estimate of u(x, t) based on observations through time t'. Error in estimating u(x, t) by the optimal estimate. E{~(x, tlt')uT(x', t" It"')} tit'; x', tit') P(x, tit; x', tit) Error covariance matrix. Any admissible estimate of u(x, t). Arbitrary vector. .. "' where u, u, w, and n are n-vectors and Pis an n x n matrix. 3 c.,•) . To facilitate working with the optimality criterion, we introduce the following inner product and norm, Y = <n(x), y(x)> = f nT(x')y(x')dDx' D = f f nT(x')E{y 1(x')y~(x")}n(x")dDx,dDx" D D IIY 11 2 = (Y, Y) = f f nT(x )E{y(x')yT(x )}n(x")dDx,dDx" 1 11 D D where y is an n-vector. We seek an estimate at time t which is unbiased and the optimal linear combination of the observations through time t. of the estimate for the case The two forms u0 = 0, h = 0, hb = 0 are w(x, tjt) = ft B(x, t, t')z(t')dt' (2.1) 0 for continuous time observations and k L B(x, tk\' tj)z(tj) w(x, tkjtk) = j =1 for discrete time observattons. the optimal B. Thus, tol fin~ the optimal w we must find The optimal B(x, t, t') Js ~enoted A(x, t, t'), i . e. G(x, t) = ft A(x, t, t 1 )z(t')dt 1 • 0 (2.2) 33 . The criterion for optimality is that the error covariance matrix be minimized. for all n. That is, we seek to minimize the scalar function Or, a~ expressed in the inner product notation; we seek " U, such that ..., II u II ~ II u - s-211 (2.3) for a 11 SG. 2.3 Summary of Results The filters derived in this chapter are summarized in Table 2.1. For all three of the observation processes, the initial and boundary conditions are " ( x , O) = u ( x ) u 0 P(x, x', O) (2.4) = P (x, 0 x') and xe:aD ( 2. 5) y + h(x, t) j =1 ~ T •GT(x', t) 1=1 -t t P(x, j=l J X., t)Mj(t)R1i(t)Mi(t)P(xi,x', t) + G(x, t)Q(x, x' t) k Instantaneous t = t Up-Date Equations - t P(x, xj' t)MJ{t)RI;(t,[z;(t) - M;(t)O(x;, t)l i =1 j at (x, x•, t) = LxP(x, x', t) + D(x, x', t)Lx' aP dA Wx, t) = LxO(x, t) + '--' Filter i =1 t T t P(x, xj, tkltk _ l)Mj(tk)Kji(tk) t)GT(x', t) • [z;(tk) - M;(tk)Q(x;, tkltk _ 1~ · j =1 O(x, tk) = O(x, tkltk _ 1) + f + G(x, t)Q(x, x', ~x, x', ti\ _ 1) = L/(x, x', tltk _ 1) + P(x, x', ti tk _ 1)L:, Predktion-Tguations ~flrx tit ) = L u(x, tltk 1) + h(x, t) . W' k-1 x tk - 1 ::s t ::s tk ~.. Discrete Time Continuous Time Observation Process Linear Filters Derived in Chapter 2 - Detenninistic Boundary Conditions and Point-Wise Observations Table 2.1 .i:,. V,J • i =1 T t L Lp P(x, xj, tk!tk _ 1)Mj(tk)Kj;(tk) p Filter j =1 [z I: • P(x;, t"ltk _ 1; x', tkltk _ 1)dt 11 dt' i=1 I: j =1 p p • P(x, tkltk _ 1; xj, t' ltk _ 1)Mj(t')K;;(tk)M;(t 11 ) ftk ftk tk - 1 tk - 1 k- 1 M1(t")a(x 1, t"I tk _ )dt"] dt' 1 i =1 L LP(x, tkltk _ 1; xj, t' ltk _ 1) 1(tk) - I:k P(x, x', tk) = P(x, x', tkltk _ 1) - • Mj(t')Kji (tk) f p T KJ1 .. (tk) = M.(tk)P(x., x., tk!tk - 1)M.(tk) + R.. (tk) J J 1 1 J1 • M;(tk)P(x;, x', tkltk _ 1) P(x, x', tk) = P(x, x', tk!tk _ 1) - .f - . t p Time -"Av~ra~e..-- " u(x, tk) = u(x, tkltk _ 1) + k ', tk - 1 j=l O,bservation Process Table 2.1 (continued) CJ, C,J Ntj' is defined as follows T • P(xj, t' ltk _ 1; xi' t"ltk _ 1)M~(t 11 )dt 1 dt 11 + Rji(tk) Mj(t') tk - 1 tk - 1 rk Filter = LxP(x, k= 1 L P 1 J lJ N. kNk* . = I 6.. tk - 1 :'.S t I < t :'.S tk The relationship between the matrices, Nik and t!tk _ 1; x', t' ltk _ 1) ~~ (x, tltk _ 1; x', t'ltk _ 1) = ~~ (x', t'ltk _ 1;x, tltk _ 1) Jl K .• ( tk) = rk The notation Rrj ~nd Krj is used in this chapter. Observation Process Table 2.1 (continued) C, vJ 37 . For the discrete-average observations, a boundary condition is also required for P(x, tltk _ 1; x', t']tk _ 1). 2.4 Derivation It is of the Filters The approach taken in this chapter to deriving the filters is based on the orthogonal projection lemma. The orthogonal projection lemma "' of the state, U, out of the set characterizes the optimal estimate, U, of admissible estimates n. The lemma states "' IIU - UII ~ IIU - nll if and only if "' ( 2. 7) (U - U, n) = 0 for all n, and furthermore, if another estimate "'U' also satisfies (2.7), then "' "' II U U' II = 0 (2.8) 2.4.1 Continuous Time Observations For the continuous time observation :Pr<2cess the optimality con1I dttion, Equation (2. 7), becomes ' ' :38 . f f 1/(x") D D r E{~(x", t)zT(t')} O (2.9) • BT(x"', t, t')dt' n(x"' )dD XudD X"' = 0 for all n and B. The following theorem leads directly to the criterion for A(x, t, t'), the optimal B(x, t, t'). Theorem 2.1 Q(x, t) satisfies the optimality condition, Equation (2.9), if and only if ' t > t' (2.10) Let B(x, t, t') = E{u(x, t)zT(t')}. Then (2.9) becomes Proof Sufficiency: Necessity: Obvious Unless B ~ 0, 8B1 will be non-negat{ve definite and thus not satisfy this equation. Thus, (2.10) holds_. # ( .J '' We will find it expedient to write Equation (2.10) as E{u(x, t)z1(t')} - ft A(x, t, t )E{z(t")zT(t')}dt 11 0 11 = 0, t > t' (2.11) 39. and (2.12) follows directly. Equation (2.10) will be referred to as the Wiener- Hopf equation (for continuous time observations). In a sense, the filtering problem is solved; we have in Equation (2.11) an implicit equation for A. However, the fonn of Equation (2.11) is impractical in that it is not recursive. more practical, We proceed to find the recursive filter. Consider the derivativ~ of the Wiener-Hopf equation with respect to time. Using the fonn in Equation (2.11) we have -r aA(x,ai' tll) E {z(t")zT(t')} dt" = 0 t > t I• (2.13) 0 After rearrangement, use of (2.11), (1.7) and (1.1), and attention to the stochastic properties of v and w, Equation (2.13) becomes (2.14) We now show that 40 . LXA ( X, t' t I ) - A( X ' t ' t) f D M( X I ' t) A( X I ' t ' t I ) d DX I 3A(x, t, t') at = 0 • t >t' • (2.15) Proof of Equation (2.15) Let C(x, t, t') denote the left hand side of (2.15). Then if A(x, t, t') satisfies the Wiener-Hopf equation, so does A+ C. Thus both u(x, t) = ft A(x, t, t')z(t')dt' and 0 a•(x , t ) = t[ A( x , t , t • ) + c(x , t , t , )J z ( t , ) dt • 0 are optimal estimates of u(x, t). By the orthogonal projection lemma (Equation (2.8)) we have But E { z(t 11 )zT(t111 )} can be shown to be positive-definite. Thus C = 0, and (2.15) holds. # Equation (2.15) is used to find a differential expression for O(x, t). Taking the derivative with respect to time of the defining equation for 4 1. u(x, t), we have Substituting (2.15) into (2.16)yields, (2.17) The gain, A(x, t, t), remains to be determined. We will make a notational change and denote A(x, t, t) as K(x, t). Starting again with the Wiener-Hopf equation and using (2.12) and (1.7} we have (2.18) which reduces to J/{d(x, t)i/(x', t')}MT(x', t')dDx' = A(x,_t, t')R(t ( ] ,' 1 ) , t > t'. (2.19) ' Both sides of (2 . 19) are continuous in t, · thus we can take the limit t' ~ t. Doing so results in 12 . (2.20) The estimate 1s specified by Equations (2.17) and (2.20) and the initial and boundary conditions O(x, o) = 0 ,.. (2.21) f\u(x, t) = 0 , xe::ao. (2.22) We can now derive the equation for the error covariance matrix. The error, u(x, t),is given by ~ (x, t) = Lxu(x, t) - K(x, t) Jo M(x', t)u(x', t)dDx' + G(x, t)w(x, t) - K(x, t)v(t) t > 0 , xe::D . (2.23) With initial and boundary conditions Su(x,t)=O ,xE:aO. (2.24) The solution to this partial differential equation can be expressed in tenns of the Green's function matrix~. as follows (2.25) 43 . where - K(x, t) f0 M(x", t)<1>{x t; x', t')dDx" 11 , limt + t•+<I>(x, t; x', t') = Io(x - x') (2.26) <I>(x, t; x', t') = O , xE:aD . Directly from the definition of P(x, x', t) and (2.26) we can determine P to be given by - K(x, t)R(t)KT(x', t) + G(x, t)Q(x, x', t)GT(x', t) (2.27) (2.28) P{x, x', t) = 0 , )(E:ao ( 2. 29) . The complete filter is given by Equations (2 . 17), (2.21), (2.22), (2.20), {2.27), (2.28), and (2.29). These i equations can be rewritten in terms of the individual components of the observation matrix and Ras follows (the initial and boundary conditions remain unchanged), 44 . au (x, t) = Lxu(x, t) + ff f0 P(x, x 11 , =1 i =1 j 111 111 (2.30) • {zi (t) - JD Mi (x , t)u(x , t)dDX111 }dox" ~~ (x, x', t) = LxP(x, x', t) + P(x, x t)LJ, 1 , p 1 + G(x, t)Q(x, x', t)G (x', t) - p J0J0 P(x, x", t ~ [ j =1 • M} (x 11 , j =1 t)Rt (xj, xi, t)Mi (x 11 : t)P(x 111 , x', t)dDxu~Dx"· (2.31) 2.4.. 2 Discrete-Instantaneous Time Observations * In considering discrete time observations we divide the problem into two parts, prediction and up-dating. During the interval tk _ 1 < t < tk no new infonnation is available, thus we pred-ict an estimate of the state based on the past estimate at time tk _ 1. At time tk an observation is recorded and we up-date the estimate based on this observation. The prediction equations were derived by Tzafestas and Nightingale (1968a) and are given in Table 2.1 . We derive the up-date equations. Using difference equations instead ~f differential equations, the I procedure of Section 2.4.1 for continuous time observations is followed here for discrete-instantaneous observa~ions. From the orthogonal pro- *The orthogonal projection lemma is applied in Appendix B to discrete time observations of lumped parameter systems. 15 . jection lerrma, the optimality condition, Equation (2.7), becomes k ff r/(x )E{u(x",tk) zT(ti)}BT(x"',tk~ t;)n(x"')dDx dDx"' = 0 , ( 2. 32) DD L 11 ) 11 i = 1 for all n and B. The theorem for the Wiener-Hopf equation is now stated . Theorem 2.2 ~(x, tk) satisfies the optimality condition, Equation (2 .32), if and only if (2.33) Proof Follows directly from Theorem 2.1. # Again, Equation (2.33) will be referred to as the Wiener-Hopf equation (for discrete-instantaneous observations), and the following forms will also be useful. L k 1 E \u(x, tk)z (t;)l - 1 A(x, tk' tj )E (z(tj)z (t;)l = 0 , j = 1 1~ i ~ k ( 2 . 34) 46 . Consider the difference of the Wiener-Hopf equation at times tk and tk - l' Using (2.34) and rearranging terms we have E { [u(x, tk) - u(x, tk _ 1)]zT(ti)} - A(x, tk' tk)E { z(tk)zT(ti)} -L k- 1 j [A(x, tk' tj) -A{x, tk- l ' tj)] E/z{tj)/{t 1 )j =1 = 0 , 1~ i ~ k- 1 . ( 2. 37) The state, u(x, tk),can be expressed in tenns of the Green's function matrix, 1/J, as follows (2.38) where : \ (x , t ; x ' , t ' ) = L i ( x_, t ; x 1 , t 1 ) ! '' 1+1/J(x,t; x', t') = Io(x · - ..x') f I lim t + t 1/J ( x, t; x 1 , t 1 ) = 0 , x Ed0 . (2.39) 47 . After rearrangement, use of (2.38), (2.34), and (1.8), and attention to the stochastic properties of v and w, Equation (2.37) becomes k- 1 L {[ft(x, tk; x', tk _ 1)A(x', tk _ 1 , tj)dDx' - A(x, tk _ 1 , tj)] j =1 - A(x, tk, tk) f0f0M(x', tk)iJl(x', tk; x 11 , tk _ 1)A(x 11 , tk _ l' tj)dDxaadl\a In a proof analogous to that for Equation (2.15) we can use the orthogonal projection lemma ( Equati on (2.8)) to prove the following (2.40) Taking the difference between G(x, tk) and u(x , tk _ 1) and sub, stituting (2.40), we find a recursive equ~tion for u(x, tk), 48 . (2.41) where O(x, tk !tk _ 1) is the previously defined optimal prediction of u(x, tk). A(x, tk' tk) is determined by substituting Equations (2.41), (2.35), and (1.2) into the Wiener-Hopf equation for i = k. Renaming A(x, tk' tk) to be K(x, tk) we obtain ·[1 M(x", tk}[u(x", tk} - il(x", tkltk - 1l]dox" + v(tk}r} = 0 (2.42) and thus, K(x, tk) : r P(x, x 111 11 ', tk ltk _ 1)MT(x , tk)dOX111 JD •(foJ0 M(x', tk}P(x', x", tkl~k _ 1}MT(x", tk} • dD X dO X I II + R( tk),·} . -1 . , ] , (2.43) The up-dated error covariance matrix is determined by direct substitution. It is 49 . (2.44) The up-dated filter is given by Equations (2.41), (2.43), and (2.44). These equations can be rewritten in terms of the individual components of the observation matrix and Ras follows, p p x'", tkltk _ 1) [: [: MJ(x"', tk)K;i(tk) j=li=l - J0M1(x", tk)il(x", tkltk _ 1)dDx" jdox'" p (2.45) p - fofoP( x, x"', tk Itk - 1) [ : [: M} (x 111 , tk) j =1 i = 1 (2.46) f Kji = J Mj(x", tk)P(x", x"', tkltk _ 1)Mr~x"', tk)dDx"dOx"' + Rji(tk) • 00 I (2. 47) 2.4.3 Discrete-Average Time Observations For the discrete-average time case the filter derivation proceeds 50 . essentially as in Section 2.4.2. The differences that do exist arise because z(tk) contains information from the entire interval tk _ 1 t ~ tk that has been compressed into one datum. ~ The derivation will be . sketched here . The Wiener-Hopf equation for discrete-average time observations is identical to the one for discrete-instantaneous time observations, Equation (2.33). Equations (2.34) and (2.35) also apply. After taking the difference of the Wiener-Hopf equation at times tk and tk _ 1 we can prove the analog to Equation (2.40), which is, Q0,(x, tk; x', tk _ 1)A(x', tk _ 1, tj)dDx' - A(x, tk _ 1, tj)] - A(x, tk' tk)ftk fofo M(x', t')•(x', t'; x", tk _ 1)A(x", tk _ l' tj) tk - 1 • dDxudDx,dt' - [A(x, tk, tj) - A(x, tk _ l' tj)] = 0, 1 ~ j ~ k - 1. (2.48) Taking the difference between u(x, tk) and u(x, tk _ 1 ), after considerable rearrangement and introduction df u(x, t'ltk _ 1) we have 1 )d0x,dt') . (2.49) /j J . Again, renaming A(x, tk, tk) to be K(x, tk), we employ the Wiener-Hopf equation at i ~ k and detennine K(x, tk), f f M(x', t')P(x', t' ltk _ -1 0 0 • MT(x", t")dDx,.d(ydt"dt' + R(tk) i-l . (2.50) The up-dated error covariance matrix is (2.51) In the equations for K(x, tk) and P(x, x', tk) we need the value of P(x, tltk _ 1; x', t'ltk _ 1) for tk _1st'< ts tk. forward to determine from the prediction equation. This is straight- P(x, tltk _ 1;x 1 ,t• I tk _ 1) is given by (2 . 52) 52 . ; x 1 , t 1 ltk _ ) given by the 1 1 prediction covariance equation and boundary condition with initial condition P(x, t• ltk _ The up-dated filter can be rewritten in tenns of the i ndi vi dual components of the observation matrix and Ras follows, Ill 1 ;x ' t"'I t · k - 1 ) (2.53) t"')K~. ( tk) Jl t t J (2.54) \ K .. ( t ) = f k f f k M. ( x 11 , t 11 ) P( x~•, t II It • x , t 111 I t k-1 ) Jl k t Dt DJ k-1' k - 1 k - 1 •, 111 (2.55) [JJ . CHAPTER 3 Distributed Parameter Filters: Stochastic Boundary Conditions 54. In this chapter we present the filters for linear distributed parameter systems with stochastic boundary conditions. Such boundary-conditions contain both detenninistic and stochastic tenns and are of the fonn where hb and wb are detenninistic and stochastic functions, respectively. 3.1 Introduction A specific class of systems from those given in Section 1.3 is considered in this chapter. In most previous work, consideration of more general cases has necessitated leaving unspecified functions in the final result. Consideration of a specific class of systems enables us to derive the filters completely such that no additional functions need to be detennined for implementation. Domain, D The class of systems we .consider is - one dimensional Spatial Operator, Lx - second order partial differential operator Boundary Operator, Bx - mixed or Dirichlet. Many engineering systems are of this clas ~ and the extension to different domains, in particular, higher dimensions, ~involves no new concepts .' relative to filtering theory. i •~ ( '1 '- We shall see that the effect of the stochastic boundary conditions on the filter is independent of the observation process. (There is one interesting exception; the case of Dirichlet 55 . boundary conditions with observations on the boundary. cussed in Section 3.3.2.) This is dis- Although we do not discuss continuous time, continuous space observation processes, such processes are, as noted in Chapter 1, frequently considered in the literature. We note that the results presented in this chapter can be applied directly to such processes. Section 3.2 contains the specifications of D, Lx, and Sx and the summary of the results. Sections 3.3.1 and 3.3.2 contain the derivations of the filters for mixed boundary conditions and Dirichlet boundary conditions, respectively. 3.2 Summary of Results As discussed in Section 3.1, the system we consider is as follows au - Lxu = G{x, t)w{x, t) + h(x, t) Lu= af XED, 0< t { 3.1) where {3 . 2) and Dis the nonnalized interval O < x < 1, with initial condition u{x, 0) = u0(x), and boundary conditions 5 G. The stochastic boundary inputs, w0 and w1 , are independent, zero mean, and white with covariance matrices Q0(t) and Q1(t), respectively. coefficients a,¢, and e are nxn matrices. The The boundary conditions are considered in three cases, Case I - mixed (flux and state tenns); ¢0 and ¢1 are nonsingular; observations in the domain, D, and/or on the boundary Case II - Dirichlet (state tenn); e0 and e1 are nonsingular, ¢0 = ¢1 = O;observations restricted to the domain Case I II - Dirichlet (state tenn); e0 and e1 are nonsingular, ¢0 = ¢1 = O;observations on the boundary and in the domain. It will be clear that the filters for combinations of the above cases (e.g. ¢0 = 0, ¢1 nonsingular) follow directly from the three cases considered. For Case I, the three observation processes discussed in Chapter 2 are all considered. The partial differential equations for u(x, t) and P(x, x', t) are not affected by the stochastic conditions for u(x, t). boundary These partial differential equations are those presented in Table 2.1 for each of the respective observation processes. The single set of boundary conditions in Table 3.1 for Case I applies to all of the observation processes considered. ' Similarly, for Case II, the partial d{fferential equations for u(x, t) and P(x, x', t) are those presented •in Table 2.1 for each of the I ./ _', respective observation processes, and the sing1e set of boundary conditions for all three observation processes is given in Table 3. 1. ~ -~ hit) ' + Q0(t) 8/(t)M~\r\t) [z~(t) - M~(t)e;j 1(t)h0(t)] 8/(1, x', t) = Q1(t)[a/1, t)0'/(t)]To'(x' - 1) 80P(O, x', t) = -Q 0(t) [a 2(o, t)0 1(t)]To'(x') 1 0 lT l , . 1 1 = 0, 1 z~, M~, and R~ are defined in Section 3.3.2. b b b -1 -lT bT Ki(t) = Ri + Mi(t)ei (t)Q;(t)ei (t)Mi (t) 1 1 bT b-l [ b b 1 (t)Ml (t)Kl (t) zl(t) - Ml(t)ei (t)hl(t½ T T -1 l T 8/0, x', t) = L0 1(t) - o1(t)e 1 (t)M~ (t)K~ (t)M~(t)e 1(t)Q 1(t)Je 1 (t)a~(l, t)o'(x' - 1) I, 81U(l, t) = hl(t) + Ql(t)ei X =1 X = 0 X = 1 X = 0 X= 1 81P(l, x', t) = Q1(t)¢ 1 (t)6(x' - 1) T X = 0 0 T 80P(O, x', t) = -Q 0(t)¢ 1 (t)o(x') I, • T T -1 ] 1T 1 8aP(O, x', t) =-Loa(t) - Oa(t)00 1 (t)M~ (t)'(~ (t)M~(t)e; (t)Qa(t) e; (t)a1(o, t)o'(x') S01Jio, Case III 8 u(l, t) = h (t) 1 1 8 G(o, 0 t) = h0(t) = h1(t) 8 u(l, t) 1 Case II = h0 (t) 8 u(O, t) 0 Case I Linear Filter Boundary Conditions for Stochastic System Boundary Conditions Table 3.1 ~ ::.-.. 58 . For Case III, the results indicate that if the system boundary conditions are of the Dirichlet type then only for continuous time observations is it useful to take observations at the boundary. only be concerned with this type of observation process. Thus, we will The boundary conditions in Table 3.1 for Case III apply only for continuous time observations. The partial differential equations for O(x, t) and P(x, x', t) are those presented in Table 2.1 for this process. Discrete time boundary observations need not be considered for reasons discussed in Section 3.3.2. 3.3 Derivation of Filters The central concept in these derivations is the formal use of the Green's function form of the solution for a partial differential equation to show that a boundary inhomogeneity may be represented in the differential equation, (Stakgold, 1967; Sakawa, 1972). In the derivations, only the x = 1 boundary condition will be taken to be stochastic . The x = 0 boundary will be taken to be deterministic. The effect of the stochastic term at x = O will be obvious and is included in Table 3.1. At both boundaries we set the deterministic term equal to 0 to avoid obscurinQ the derivations . Before treating Cases I and II separately and in detail, we develop here the conmen elements of these derivations. Much of the background mathematical development is identical in t~ese derivations and they follow the same logical arguments . (The d~rivation in Case III also draws on these concepts, but remains substdnti~lly different.) In Cases I and II the derivations proceed as follows. 59 . Step 1 The system boundary conditions are made deterministic by specifying the stochastic inhomogeneity in the partial differential equation. Step 2 Application of the known filter for the case of deterministic boundary conditions yields an estimate equation with deterministic boundary conditions. The effect of the stochastic inhomogeneity enters in the covariance equation. Step 3 The effects of the stochastic inhomogeneity is specified in the boundary conditions for the covariance. As background material, we discuss first the form of the filters previously derived for deterministic boundary conditions and then the Green's function form of the solution to a partial differential equation. We start by noting that the system and filters presented in Chapters 1 and 2 may be expressed as follows, (3.5) - System - Filters - Continuous time observations au af - L ,... = S,,.(x , 't ) + H( X, XU u . I t) (3.6) j , ," ap - L p - PLT af X x' = S-p ( X' XI' t) + Ow ( x , x • , t) ( 3. 7) 60 . - Discrete time observations - Prediction (3.8) (3.9) - Up-Date Not a function of W, with appropriate initial conditions. The boundary conditions are deterministic and in particular those for the covariance equation are homogeneous. H represents all detenninistic inputs and Wall stochastic inputs, with E { W(x, t)WT(x', t')} = Qw(x, x', t)o(t - t') . s0 and Sp, are not functions of W. Equations (3.6) - (3.9) show the influence of stochasfic inputs ·on the filtering equations. Their effects are the same for all of the observation processes considered. We see that Wenters into the covar- iance equation and that the estimate equation is unaffected by the presence of the stochastic inputs. In the filter derivations we will . 'add ' stochastic inputs to the system equation, thus it is important that we see clearly what the effects will \ be. We now develop the Green ' s function form of the solution to systems of the forms Lu = :~ - Lxu = Tu( x, t) and (3.10) 61 . (3.11) where Lx and the domain are as given in Section 3.2. * defined by the Introduce the adjoint Green's function matrix, l/J, following system L: l/J*(x, t; x', t') = -~i* - L:l/J* = 0 (3.12) (3.13) with end condition lim t + t•- l/; * (x, t; x', t') = Io(x - x') (3.14) and boundary conditions to be specified later. The, operators Lx* and Lx* are the adjoints of the operators Lx and Lx, respectively. Note that l/; * is related to the Green's function matrix l/J, introduced in Equation *T (2.39), by the relationship l/; (x', t'; x, t) = l/;(x, t; x', t'). In spite of this simple relationship we prefer to use l/; * in this discussion because it is most natural to do so. By integration by parts it can be shawn that l (3.15) and Ii * r] Pl/! * - l/J.rPllx"'l/J * - [ l/J *·r Lx" *l}J dx"'dx"dt' = f 1 f 1 l/J*rPl/!*dx"'dx t 0 0 0 11 (3.16) Using Equations (3.12) and (3.14) we find the 'solution' for u and P to be u(x, t) = + l f0 l/J *T (x', O; x, t)u (x')dx' 0 t fl l/J *T (x', t'; x, t)L ,u(x', t')dx'dt' 0 0 X f - f:[-w•T(x', t';x, t)a (x', t') au(~~: t') 2 + ( ~x' (w•T(x', t'; x, t)a2 (x', t')) - iµ*T(x', t'; x, t)a (x', t•~) u(x', t')] dt' : 1 (3.17) j and P(x, x', t) = l fl l/J•T (x", O; x, t)P (x", x"' )l/!* (x 11 ', O; x', t)dx"'dx" 0 f0 0 63. JJ t 1 1 *r [ aP( x 'a~•' t, -) - Lx P{x + O O O "l/1 (x", t'; x, t) f 11 '" 11 t f1 [ .r fo o -"l/1 (x II , t, ; x, t) a2 ( x" , t, ) 11 , x"', t') clP(x", x"', t,, - ax ...L II - "l/J*T(x 11 , t'; x, t)a (x 11 , t')) P(x", x"', t')]"l/J*(x 11 ' , t'; x', t)dx 11 dt 1 l 1 , 0 - fat fal "l/1*T{x ( 'a~m' t' ) a2(x T 111 , t')"l/1* (x"', t'; x', t) t'' ; x, t) [ _aP x II 11 , "' + P(x", x"', t 1) ( ~x"' ( a~(x'", t' )"l/1* (x"', t'; x', t)) (3.18) Equations (3.17) and (3.18) hold for the boundary conditions of both Case I and Case II. The boundary conditions on "l/1* have yet to be specified. To proceed with our discussion we must now consider the two cases individually. 3.3.1 Mixed Boundary Conditions 'I For Case I, the adjoint boundary conditions are U4 . * *T BOVJ = *T *T -1 ] [ ~x (VJ a2) - VJ (al - <Po 0 o) X * *T Bl VJ = = 0 = 0 .r a *T [ 'ax (VJ a2) - VJ (al - q,llel)] (3.19) = 0 X = 1 Recall that we have set only the x = 1 boundary condition on u to be stochastic. Thus the system boundary conditions are (3.20) Combining Equations (3.17), (3.19) and (3.20) we have for u, u(x, t) = 1 *T J0 VJ (x', O; x, t Jl *T + f0 0 + t *T l O VJ (1, t'; x, t)q, (t')w1(t')dt' VJ (x', t'; x, t)Tu(x', t')dx'dt' f 1 which is equivalent to u( x , t) = l f 0 •T VJ ( x' , 0 ; x , t) u0 ( x ' ) dx ' ; . (3.21) 65 . (3.22) Equation (3.22) is the solution to (3.23) with deterministic boundary conditions. Thus we have shown that u can be equivalently described by (3.10) with (3.20) as the boundary condition or by (3.23) with deterministic boundary conditions. Based on Equations (3.18) and (3.19), we can also show that the covariance, P, described by Equation (3.11) with boundary conditions s0P(O, x • , t) = o (3.24) is equivalent to the specification , T 1 (t)o(x' - 1)2 1 + 2o(x - l)<Pi (t)Q (,t_)<P 1 I 1 1 ,' with homogeneous boundary conditions at both x = 0 and x = 1. (3.25) 66 . We now have all of the background material and can proceed through the derivation. Step 1 - The system is described by Equation (3 . 1) and the stochastic , boundary conditions (3.20). An equivalent specification of the system is (3 . 26) with deterministic boundary conditions. Step 2 - The filter for this system is known from the results of Chapter 2. It is ( 3. 27) s0u(O, t) = 0 (3. 28) f3ou(l, t) = 0 aP at -. T + 2o(x - 1)¢ 11(t)Q 1(t)¢ 11 (t)o(x' - 1)2 (3 . 29) 67 . s0P(O, x', t) = O (3.30) s0P(l, x', t) = 0. Step 3 - The covariance described by Equations (3.29) and (3.30) is equivalently specified as 9-f. at _LP X - PLTx• = Sp(x, x 1 , t) + G(x, t)Q(x, x', t)GT(x', t) (3.31) BaP(O, x', t) = 0 (3.32) B1P(l, x', t) = Q1(t)¢i 1T (t)&(x' - 1) . Equations (3.27), (3.28), (3.31) and (3.32) and the appropriate initial conditions represent the filters for the stochastic boundary conditions of Case I. 3.3.2 Dirichlet Boundary Conditions In Cases II and III we consider pi/ iJhlet boundary conditions. First we will consider Case II; that is, toe case of no observations on the boundary. 68 . The system boundary conditions are (3.33) As for Case I, Equations (3.17) and (3.18) are the main elements in this derivation and we start by defining the appropriate adjoint boundary conditions. For Case II, the adjoint boundary conditions are * *T = \jJ *T So l/J 80 =0 X = 0 {3.34) * *T = \jJ *T S1 l/J 81 X =1 =0 . Combine Equations (3.17), (3.33) and (3 . 34) and note that the following relationship holds by the properties of the derivative of the delta function {3.35) 69. It thus follows that u can be described by the equation Lu= Tu(x, t) + a2(1, t)e 11(t)w1(t)2o'(x - 1) (3.36) with deterministic boundary conditions. Based on Equations (3.18) and (3.34), we can also show that the covariance, P, described by Equation (3.11) with boundary conditions (3.37) is equivalent to the specification (3.38) with homogeneous boundary conditions. The three step derivation previ- ously outlined is now followed and the results are as stated in Table 3.1. We conclude this chapter on stochastic boundary conditions . by considering continuous time observ 9tions systems with Dirichlet boundary condition~. on the boundary for As indicated in Table 3.1, we introduce a minor amount of additiona'l.1n.otation in treating this type of observation process. Let there be p + 2~·observations with the total observation vector denoted Zand defined by 70. z(t) Z(t) = M(x 1 , t)u(x', t)dx' z0b(t) v(t) J10 b M0( t )u( 0, t) = + M~ ( t) u(1 , t} z~(t) v~(t) (3.39) v~(t) where z represents the p observations in the domain and z~ and z~ represent the boundary observations at x = 0 and x = 1, respectively. The boundary observation errors, v~ and v~, are zero mean and have second moment properties defined by v( t) E v~(t) [vr(t ' ) v~\t') v~\t•)] R(t) 0 0 0 R~(t) 0 0 0 R~( t) v~(t) = tS(t-t') (3.40) where R(t) , R~(t), and R~(t) are positive definite, and we set v~ and v~ to be independent of each other ·and of v. The terms z, M, v and R are defined as they were in Chapter 1. In deriving the filter for this case we will only consider the x = 1 boundary to be stochastic . Thus we will t1ot consider observations at the ( ' ' ' x = 0 boundary and the terms fo Equations, (3-.39) -and (3.40) relating to x = 0 ·are to be removed. The system boundary conditions are 71. {3.41) Finally, as in the derivations in Chapter 2, we set h = 0 and u = O. 0 The system may now be specified with deterministic boundary conditions and system equation We seek the estimate of the fonn w(x, t) = J: {3.43) [B(x, t, t') which minimizes the error covariance matrix. Following the conventions in Chapter 2, u, A, and A~ are the optimal w, B and B~, respectively. Based on the orthogonal projection lemma the optimality condition is T ] [BT (XII, t, t I) ] zb ( t') T l B~ (x 11 , t, t') • dt'n(x 11 )dx 1 dx 11 = ·o I l , (3.44 ) ,j •~-, for all n, Band B~. This condition leads dj.rectly to a Wiener-Hopf equation which can then be expressed as the following two separate relationships ?2 . (3.45a) and E{u(x, t)z 1bT (t') } = 0 , t > t' . (3.45b) Alternative expressions for these equations are + A~(x, t, t )E{z~(t )zT(t 11 11 and f: E\u(x, t)z~T(t') )- 1 )}] dt 11 = 0, t > t' (3.46a) [A(x, t, t")E\z(t")z~\t•)) (3.46b) Differentiating (3.46a) and (3.46b) with respect tot we have J: C(x, t, t")E(z(t")zT(t')) dt" + J: C~(x, t, t")E(z~(t")zT(t'))dt" = 0, t > t' (3.47a) and T t f C(x, t, t )E{z(t )z~ (t')} dt" 11 11 0 (3.47b) where 73. C(x, t, t") = - aA(x,ai' t"} + LxA(x, t, t 11 ) (3.48a) - A(x, t, t)J: M(x', t)A(x', t, t")dx' and aA~{x, t, t~) b • at + LxAl ( x ' t , t" ) (3.48b) Equations (3.48a) and 3.48b) are not exactly analogous because, for t > t', E{z~{t)z~~t')} and E{z~{t)zT(t')} vanish whereas E{z(t)zT(t')} and E{z{t)z~T(t')} do not. It can be shown that C and C~ vanish and· that the differential expression for u becomes (3.49) with deterministic boundary conditions. A(x, t, t) and A~(x, t, t) are found to be A(x, t, t) = K(x, t) = J1 P(x, x 0 1 , t~M1 (x 1 , t)dx'R- 1(t) ; j (3.50) ' b T. T -1 Al(x, t, t) = a2(1, t)ei 1(t)Ql{t)ei 1 (t)M~ (t)K~ (t)2o (x - 1) 1 (3.51) 74 . where (3.52} The error covariance equation is thus, ap ( X, X , t = af I } (3.53} with homogeneous boundary conditions. The filter can alternatively be written with an inhomogeneous boundary condition at x = 1. The partial differential equations for "u and Pare then those in Table 2.1 and the boundary condition is that given in Table 3.1. The filters for discrete-instantaneous time and discrete-average time boundary observations have also been investigated. discrete time filters are impractical. Both of these In the instantaneous case no ' infonnation is gained by taking such observations. In the case of averaged observations, the computational burden would be excessive. In the instantaneous case, the reaso~ / br this result is clear. Discrete-instantaneous time boundary observations yield a characterization of an input (the boundary value) at an instant. The influence of an ? uc . input at an instant is infinitesimal; thus, no infonnation is gained that can be used to estimate the state. 7(j . CHAPTER 4 Summary 77. In this work, we have considered the linear distributed parameter filtering problem for systems of engineering interest. Fonnal proce- dures have been employed in developing filters for a variety of obser- vation processes and boundary conditions. Integral and point-wise spatial observation processes, along with continuous and discrete tenporal observation processes, have been considered. For the case of deterministic boundary conditions, the filters have been developed for a general linear system. In order to provide filters which are completely specified, for systems with stochastic boundary conditions a specific linear system has been considered. ?8. References Aidarous S. E. and Ghonaimy M.A. R. (1976) Optimal filtering for a class of stochastic distributed systems. Int. J. Systems Science, ?_, 545 - 556. Ajinkya M. 8., Ray W. H., Yu T. K. and Seinfeld J. H. (1975) The application of an approximate non-linear filter to systems governed by coupled ordinary and partial differential equations. Int. J. Systems Science,§_, 313 - 332. Athans M. (1970) Toward a practical theory for distributed parameter systems. I.E.E.E. Trans. Auto. Control, AC-15, 245 - 247. Atre S. R. (1972) A note on the Kalman-Bucy filtering for linear distributed-parameter systems. I.E.E.E. Trans. Auto. Control, AC-17, 712 - 713. Atre S. R. and Lamba S.S. (1972a) Derivation of an optimal estimator for distributed parameter systems via maximum principle. I.E.E.E. Trans. Auto. Control, AC-17, 388 - 390. Atre S. R. and Lamba S.S. (1972b) Optimal estimation in distributed processes using the innovations approach. I.E.E.E. Trans. Auto. Control, AC-17, 710 - 712. Atre S. R. and Lamba S.S. (1972c) Filtering for linear distributed parameter systems via boundary measurements. Proc. I.E.E., 119, 757 759. Balakrishnan A. V. (1974a) Stochastic optimization theory in Hilbert spaces - I. J. Appl. Math. and Opt.,!., 97 - 120. Balakrishnan A. V. (1974b) Identification and Stochastic Control of a Class of Distributed Systems with Boundary Noise. University of California, Los Angeles, Rept. UCLA-ENG-7461. Balakrishnan A. V. and Lions J. L. (1967) State estimation for infinitedimensional systems. J. Computers and Systems Science,!., 391 403. Bensoussan A. (1971) Filtrage Optimal des Systemes Lineaires, Dunod, Paris. ' Curtain R. (1975a) A survey of infinite-di,~n$ional filtering. 1, ', Review, !.Z., 395 - 411. S.I.A.M. •, Curtain R. (1975b) Infinite-dimensional filtering. g, 89 - 104. S.I.A.M. J. Control, 7.9 . Falb P. L. (1967} Infinite dimensional filtering: the Kalman-Bucy filter in Hilbert space. Inf. and Control,!!., 102 - 137. Fujishige S. (1975) Minimum-variance estimation for a linear continuousdiscrete system with noisy state-integral observation. I.E.E.E. Auto. Control, AC-20, 139 - 140. Hwang M., Seinfeld J. H. and Gavalas G. R. (1972) Optimal least-square filtering and interpolation in distributed parameter systems. J. Math. Anal. Appl., 39, 49 - 74. Jazwinski A.H. (1970) Stochastic Processes and Filtering Theory. Academic Press, New York. Kalman R. E. and Bucy R. S. (1961) New results in filtering and prediction theory. J. Basic Eng., Series D, 83, 95 - 108. Koda M. (1975) Dynamic Estimation of Diffusion Systems-Application to Prediction of Air Pollution, Ph.D. Thesis, University of Tokyo, Tokyo. Kumar G. R. V. and Sage A. P. (1972) An innovations approach to distributed parameter filtering and smoothing. Inf. Science, i, 193 - 215. Kumar S. and Seinfe.ld J. H. (1978) Optimal location of measurements for distributed parameter estimation. I.E.E.E. Trans. Auto. Control, AC-23, To appear. Kushner H.J. (1970) Filtering for linear distributed parameter systems. S.I.A.M. J. Control, 8, 346 - 359. Lamont G. B. and Kumar K. s. P. (1972) State estimation in distributed parameter systems via least squares and invariant imbedding. J. Math. Anal. Appl., 38, 588 - 606. Lee K. Y. (1976) Optimal estimation of first order distributed systems: multiplicity and martingale approach. Preprints Joint Auto. Control Conf., 634 - 636. Meditch J. S. (1970) On state estimation for distributed parameter systems. J. Franklin Inst., 290, 49 - 59. - ' Meditch J. S. (1971) Least-squares filteri ng and smoothing for linear distributed parameter systems. Automatica, ]_, 315 - 322. Padmanabhan L. and Colantuoni G. (1974) sequertial estimation in distributed systems. Int. J. Systems st,ence, _§_, 973 - 986 . Phillipson G. A. (1971) Identification of Distributed Systems, Elsevier, New York. 80. Polis M. P. and Goodson R. E. (1976) Parameter identification in distributed systems: a synthesizing overview. Proc. I.E.E.E., 64, 45 - 61. Ray W. H. (1977a) Some applications of state estimation and control theory to distributed parameter systems. Control Theory of Systems Governed b Partial Differential E uations (Edited by Aziz A. K., . W. an a as M. J . . Academic Press, New York, 231 Ray W. H. (1977b) Some recent applications of distributed parameter control theory--a survey. 2nd I.F.A.C. Symposium on Control of Distributed Parameter Systems, Coventry. Ray W. H. and Seinfeld J. H. (1975) Filtering in distributed parameter systems with moving boundaries. Automatica, .ll, 509 - 515. Reddy K. S. and Rajamani V. S. (1975a) Optimum filtering for a class of distributed-parameter systems with coloured measurement noise. Int. J. Systems Science,§_, 615 ~ 620. Reddy K. S. and Rajamani V. S. (1975b) Optimal estimation in linear distributed systems with time delays. Int. J. Systems Science,§_, 859 - 863. Robinson A. C. (1971) A survey of optimal control of distributed-parameter systems. Automatica, Z., 371 - 388. Sahgal R. and Webb R. P. (1972) Nonlinear distributed parameter estimation. Proc. 1972 I.E.E.E. Conf. Decision and Control, New Orleans, 89 - 93. Sakawa Y. (1972) Optimal filtering in linear distributed-parameter systems. Int. J. Control, 1.§_, 115 - 127. Seinfeld J. Ho, Gavalas G. R. and Hwang M. (1971) Nonlinear filtering in distributed-parameter systems. J. Dynamic Systems, Measurement, and Control, 93, 157 - 163. . Sherry H. and Shen D. W. C. (1973) Combined state and parameter estimation for distributed-parameter systems using discrete observations. Identification and S stem Parameter lEstimation (Edited by Eykhoff P. . North-Ho an, sterdam, 6 1 ~ 6 . Shukla V. and Srinath M. D. (1972) Optima, ·filtering in linear distributed-parameter systems with multiple . time delays . Int. J . Control, 1.§_, 673 - 687. Stakgold I. (1967) Boundary Value Problems of Mathematical Physics , vol. 1 and 2. Macmillan, New York. 81. Thau F. E. (1969) On optimum filtering for a class of linear distributed-parameter systems. J. Basic Eng., Series D, 21, 173 - 178. Tzafestas s. G. (1969) Boundary and volume filtering of linear distributed-parameter systems. Electronic Letters,~, 199 - 200. Tzafestas S. G. (1972a) Bayesian approach to distributed-parameter filtering and smoothing. Int. J. Control, 12_, 273 - 295. Tzafestas s. G. (1972b) On optimum distributed-parameter filtering and fixed-interval smoothing for colored noise. I . E.E . E. Trans. Auto . Control, AC-17, 448 - 458. Tzafestas s. G. (1973) On the distributed parameter least-squares state estimation theory . Int. J. Systems Science, _i, 833 - 858 . Tzafestas S. G. (1976) Nonlinear distributed-parameter filter ing using the Fokker- Planck equation approach. J. Franklin Inst. , 301 , 429 450. Tzafestas s. G. (1978) Distributed parameter state estimation . Distributed Parameter S stems - Identification , Estimation and Control Edited by ay . H. an a n ot,s . •l . arce New York, 135 - 208. Tzafestas s. G. and Nightingale J.M. (1968a) Optimal filtering , smoothing and prediction in linear distributed-parameter systems . Proc. I.E.E., 115, 1207 - 1212. Tzafestas s. G. and Nightingale J . M. (1968b) Concerning opt imal filtering theory of linear distributed-parameter systems. Proc . I.E.E., 115, 1737 - 1742. Tzafestas s. G. and Nightingale J.M. (1969) Maximum-likelihood approach to the optimal filtering of distributed-parameter systems . Proc . I. E.E. , 116, 1085 - 1093. Tzafestas s. G. and Nightingale J . . M. (1976) Stochastic distr ibut ed filtering approach to gamma-ray imaging . Int . J . Systems Science , ?_, 1249..: 1274 . ' • f e 1d J . H. and Ray W. H. (1A ~ 74) Fi ltering in nonlinear Yu T. K. , Se,n time delay systems . I . E. E. E. Trans . Auto . Contr ol . AC-19, 324 333 . APPENDIX A Orthogonal Projection Lemma 83 . In this appendix we provide a proof of the orthogonal projection lerrma. The lemma is well known and the proof is essentially that of Kalman and Bucy (1961). Because the proof is available elsewhere, we include it here only as a reference-point for the reader. This being our purpose, we do not use the abstract fonn of the lemma as it is stated in the introduction to Section 2.4. Rather, we consider it as it specifically applies to distributed parameter estimation. Lemma nT(x')E{[u(x', t) - u(x', t)] [u(x", t) - O(x", t)]T}n(x")dD ,dD" JI DD X X sf DfDnT(x')E{[u(x', t) - w(x' , t)][u(x", t) - w(x", t)]T}n(x")dD ,dD X 11 X (A.1) (i) if and (ii) only if for all n and w,and (iii) furthennore, if another estimate a2 also satisfies (A.2) then = 0 . (A. 3) 84. Proof (i) Consider the identity +2J f nT(x )E{[u(x DD 1 1 , t) - O(x', t)][O(x t) - w(x 11 , 11 , t)]T}n(x")dD ,dD . X 11 X Note that O - w is an admissible estimate of u, thus by the hypothesis, Equation (A.2), the middle tenn vanishes and (A.1) holds. (ii) Assume (A.2) does not hold. Then there exists an w1 such that Consider another estimate w2 = u + crw1. Then ' 85 . - 2ap + a 2 f f nT(x')E{[w1(x', t)] [w1(x" , tf}n(x )dDx,dDx" • 0 0 11 For an appropriate a, the sum of the last two tenns will be negative and thus (A.1) will not hold. (iii) If both Q and u2 satisfy (A.2), then Note that Q - a2 is an admissible estimate of u, thus (A.3) follows. # 86 . APPENDIX B Certain Lumped Parameter Filters 87. In the course of developing the linear distributed parameter filters presented in this work, certain lumped parameter filters were also developed. First, we consider discrete time observation processes with con- tinuous systems. The resulting filters are not new results. However, the derivations are presented here because they represent applications of the orthogonal projection lenma technique and as background material to the distributed parameter filters in Chapter 2. Second, we discuss the influence on a lumped parameter filter of a determi ni st,i c input and a nonzero mean initial condition. The results are entirely intuitive and, in fact, in the distributed parameter developments we dismiss them as such. We feel that it is of value to present, at least in a simple case, a proof of the intuitive results. B.l Introduction Consider the system governed by the linear ordinary differential equation du (t) = F(t)u(t) + G(t)w{t) + h(t) dt defined fort> 0. (B.l) The state vector, u, is n-dimensional. w(t) is an ' n-dimensional stochastic input, white in tfme, his a known n-dimensional input, and F(t) and G(t) are known nxn matrkes. for (B.l) is The initial condition ' J '' (B.2) 88 . When the stochastic tenns are absent, the problem is assumed to be well-posed. The initial condition is unknown, with its first two moments given, (B.3) The stochastic input has the properties E{w(t)l = 0 (B.4) where Q is a non-negative definite nxn matrix. The three types of observations considered are continuous z(t) = M(t)u(t) + v(t) (8.5) discrete-instantaneous (B . 6) discrete-average . z(tk) = ftk M(t')u(t')dt : + v(tk) . (B . 7) tk - 1 The observation vector, z, is m-dimensiondl ~· ,Mis the known mxn observation matrix and vis them-dimensional observation noise. ous observations the properties of v(t) are For continu- 89. E{v{t)} = 0 {B.8) where R(t) is a known positive definite mxm matrix. Similarly, for discrete time observations the properties of v(tk) are (B.9) where R(tk) is a known positive definite mxm matrix. Finally, u0 , w and v are all defined to be independent of each other. In working with the filters we will use the following notation Optimal estimate of u(t) based on observations through time t'. O(tlt') O(t) = O{tlt) u(tlt') = u(t) - O(tjt') P{tlt'; t" It"~ = E{a< t It ' ) P(tlt') = P{tlt'; tit') P{t) = P{tlt; tit) w(tlt') Error in estimating u{t) by the optimal estimate. arc t " 1t "~,} I ) Error covariance matrix. ~Any admissible estimate of u(t). 90. n Arbitrary vector. y -- nTy = E{<n, y ><n, y >} . -_ nTE{Y1Y2 T} n 1 2 = (Y, Y) = nTE{yyT}n = <n,y> (Yl' Yz) 11 v11 2 where 0, u, w, n and y are n-vectors and Pis an nxn matrix. The fonn of the estimate we seek will be defined in each of the following sections . The criterion for optimality is that the error covariance matrix be minimized. That is, we seek to minimize the scalar function for all n . Or, as expressed in the inner product notation, we " seek U, such that (B.10) for all n. The optimal estimate, "U, is characterized by the orthogonal projection lelllTla, which is given in Section 2.4. In Sections B.2 and B.3 we consider di;screte-instantaneous and discrete-average observations, respectively. As in the distributed param- eter cases we derive these filters under , ~he,' conditions h = 0 and uO = O. Then, in S.ection 8.4 we discuss the effects of removing these restrictions in the case of continuous observations. 91. B.2 Discrete-Instantaneous Observations This problem has been discussed by others, for example, see Jazwinski (1970). Of course, this problem must be considered in two parts; prediction in between observations, and up-date at the time of an observation . The prediction equations were originally der ived by Kalman and Bucy (1961) and are given below. We consider the derivation of the up-date equations. The prediction equations are (8.11) (8.12) The up-dated estimate which we seek is of the fonn k w(tkltk) = L (B.13) B(tk' tj)z(tj) j =I Thus to find the optimal w we must find the optimal B. denoted A, i.e. The optimal Bis k L A(tk, tj)z(tj) . j =1 Using difference equations instead of ·di fferential equations, the procedure of Kalman and Bucy (1961) for continuous observations is followed here for discrete-instantaneous observations. From the orthogonal projection lerrma, the optimality condition, Equation (2.7), becomes k nT L E{u(tk)zT(tk)}BT(tk' ti)n 0 , (8.14) = i =1 for all n and B. The Wiener-Hopf equation (for discrete observations) follow directly. The estimate must obey the conditions (B.15) Now consider the difference of the Wiener-Hopf equation at times tk and tk - 1 (B.16) The state, u(tk), can be expressed intenns of the system transition matrix,$, as follows u(tk) = w(tk;tk _ l)u(tk _ 1) + Itk w(tk; t 1 )G(t 1 )w(t 1 )dt 1 (B . 17) tk - · 1 where ~ (t; t') = F(t)w(t; t') (B.18) $( t; t) = I . 93 . Using this expression for u(tk), Equation (B.16) becomes k- 1 C(tk; tj)E{z(tj)zT(t;)} = 0 , 1 5 i 5 k - 1 j =1 L where (B.19) (B.20) C(tk; tj) vanishes and we are then able to find a recursive equation for Q( tk) ' (B.21) where O(tkltk _ 1) is the previously defined optimal prediction of u(tk). A(tk, tk) is detennined by using the Wiener-Hopf equation for i = k. The up-dated error covariance matrix follows by direct substitution. Together with (8.21), the following equations comprise the up-dated filter A(tk' tk) = K(tk) . , Itk _ )MT(tk) + R(tk)j,l-1 = P(tk Itk _ 1)MT(tk) [ M(tk)P(tk 1 (B.22) 94. (B . 23) B.3 Discrete-Average Observations Using the innovations technique, this filter has also been derived by Fujishige (1975) . Because the derivation using the orthogonal projection lemma so closely parallels derivations already presented, we will make only a few co1T1T1ents here. The appropriate Wiener-Hopf equation is (B.15). From this it follows that the difference equation for A is given by [ ~(tk;tk _ 1)A(tk _ !' tj) - A(tk _ !' tj~ - A(tk, tk) ftk M(t')~(t';tk _ 1)dt'A(tk _ l' tj) tk - 1 - [ A( tk, tj) - A( tk _ l, tj)] • 0 , P j < k - !. (B.24) The filter equations are O(tk) = O(tk,tk _ 1) + K(tk)[~(tk) _ rk M(t')O(t• ,tk _ l)dt•J tk - 1 (B . 25) (B . 26) 95 . ( B. 27) along with which is determined from the prediction equations. B.4 Continuous Observat-ions wlth .a .Determin·isti'c 'Input and N-on-Zero Mean Initial Condition In this section we discuss the influence on a lumped parameter filter of a deterministic input and non-zero mean initial condition . We presen t a proof of the intuitive results. Of course, the results are that the initial condition for the estimate is the mean of the system initial con dition and that the deterministic input ap~ears as an additive term in the differential equation for the estimate . . The essential distinction between the ~e~tvations previous ly presented and this one is the form of the estimate . Instead of seeki ng an estimate which is only a linear combination of the observations ; the 96 . estimate in this case also incorporates our a priori knowledge of the influence of hand the form O(t) • r: u0 on the state. Thus, the optimal estimate is of A(t, t')z(t')dt' + r: A1(t, t')h(t')dt' + A0(t)U0 • (B . 29) In addition to requiring that the estimate minimize the error covariance, we also require that it be unbiased {i.e., E{O(t)f = E{u{t)l). Imposing this requirement constrains A1 and A0 , and the estimate takes the fonn 0(t) • r: -r:. A(t, t')z(t')dt' + r: [ $(t; t') A(t, t")M(t")$(t"; t')dt"] h(t')dt' + [$(t; 0) - r: A(t, t' )M(t' )$(t'; O)dt,J u0 where~ is the system transition matrix. {8.30) Note that the estimation error can be expressed as ii(t) • u(t) - 0(t) • [$(t; 0) (u(0) - -r: t' + fo "ol + r: $(t; t')G(t')w(t')dt~ A(t; t•)[M(t•)&(t'; O)(u(p) - u0 ) ·, 7] dt' . ~(t'; t 11 )G(t 11 )w(t 11 )dt 11 t ,:V (t')j Now consider the state u* defined by {8.31) .9 7 . (B . 32) where the initial condition has the properties (B . 33) and F, G, w, and P0 are as previously defined. Consider also the observation (B . 34) and the estimate (B .35) u* (t) -- Jto A*(t , t ' )z*(t ' )dt ' . A Note that u*(t) = ijJ(t ; O)(u(O) - Uo) + ft ijJ(t, ,t ' )G( t' )w(t ' )dt ' , JQ (B.36) ' i and that the estimation error for u* is the .~ame as that for u, wi th A 1 ' replaced by A*. Thus , A(t , t') = A*(t , t 1) , and it fo l lows tha t ~ O(t) = o.(tl + J: iµ(t, t') h(t ' )dt ' + w(t, oJ~ 0 • (B . 37) 98 . We know O*(t), thus the filter follows directly, i (t) F(t)O(t) + K(t)[z(t) - M(t)O(t)] = + h(t) (B.38) Q(O) = u0 (B.39) p(O) = Po K(t) = P(t)MT(t)R- 1(t) . 99. PART II. STATISTICAL ANALYSIS OF AIR POLLUTANT OBSERVATIONS AND MODEL PREDICTIONS 100 . Chapter 1 Introduction 101. Efforts to analyze urban air quality are often beset by the seemingly contradictory problems of a simultaneous paucity and abundance of pertinent inf?rmation. The required meteorological, emissions and ambient concentration data reach virtually unmanageable proportions. But more data, of even finer spatial and temporal resolution, are still needed to answer unsolved questions. To deal effectively with both of these problems, techniques that are essentially statistical must be introduced. Statistical techniques allow large volumes of information to be summarized compactly by identifying basic characteristics of the data. Then, knowledge of these basic characteristics allows the use of relatively small amounts of further information to interpret future or related data sets. Observations and model predictions serve complementary purposes in engineering. The observations form the ultimate base on which models are built, and the models help in understanding the principles that lead to the observations. In Chapter 2 the role of observations in assessing model accuracy is discussed. Here, great amounts of information have to be summarized by relatively small number of measures of performance. The need is to be able to make a decision as to the adequacy of a given model. In Chapter 3 observed frequency distributions are investigated through the use of highly simplified models. Here, the underlying ~ statistical nature of air quality concentrations is explored. 10 2 . Chapter 2 Evaluation of Air Quality Models: The Accuracy of Predictions Relative to Observations 103. 1. INTRODUCTION The development and application of air quality simulation models ar e areas of fundamental interest in the effort to control air quality in urban airsheds. Air quality modeling has reached the stage of maturity at which a large variety of increasingly sophisticated models are being developed and proposed for applications . Thus there exists a need to concurrently develop a framework for the evaluation of these models . Then, within such a framework, consistent and general schemes need to be developed to perform the various model evaluation tasks. Model evaluation is a problem of long-standing concern, and nearly every discussion of a given model's development or application contains comments relative to some aspect of that model's performance. In Table 1 representative references to prior air quality model evaluation studies are listed. A clear need exists for a comprehensive discussion of air quality model evaluation along with the development of general evaluation methods . The research presented here is directed at two aspects of this probl& First, in the remainder of this section, we develop a fr~mework for the eva 1uation of the performance of air quality models. Second, in Sec ti on 2 we present quantitative methods for assessing the accuracy of model pre' dictions relative to physical obs~rvations. • The final result is a package of accuracy assessment methods thit when exercised will provide a meanin ' ful evaluation of this aspect of 1 tjhe_, performance of an air quality model. This package , AQMAAP , has been coded in FORTRAN. code is given in the Appendix . A discuss ion of the Comparison of three models. Measures of correlation, RMS error, and goodness-of-fit. Plots of scatter and residuals. Evaluation of a single model. Graphic presentations of time series, frequency distribution, correlation and bias. Comparison of three models. Regression analysis. Liu et al. (1976) Anderson et al . (1977) Maldonado and Bullin (1977} Measures of bias and spread. Comparison of three models. Measures of correlation and RMS error. Regression analysis. - ....~- Johns~n et al. (1976) ,, Statistically oriented discussion. Measures of bias, spread and correlation. Regression analysis. Brier (1975) ..-:.....:,._ Evaluations Discussed Comparison of nine models. Measures of correlation and bias (based on ratios). Reference Nappo (1974) Air Quality Model Evaluation Table 1 .c,. N C) 105 . 1.1 Evaluation Framework The evaluation of an air quality model generally focuses on the accuracy of the model s predictions relative to observations. 1 This is the approach taken in this work, and it is a logical approach because it attempts to answer the 'bottom line' question, 11 predict actual concentrations? or Does it work? 11 11 How well does the model 11 However, it is often overlooked that this approach represents only a part of the total model evaluation problem . Two other major issues must be addressed: the first is validity and the second is efficiency . Before focusing on model accuracy in the next section, we will propose an overall evaluation framework, the intent of which is to provide a perspective for accuracy assessment. In Figure 1 we outli.ne the model evaluation problem from the stage of fundamental model development through to application. Validity refers to the fundamental correctness of the model formulation . The assessment of model validity is based on scientific pri ncipl es, t.e. the relationships or equatinns which comprise the model . The fi r st issue is the adequacy of the r epresentation of the physics and the chemi stry of the system . When the system is complex, it is difficult to index quantitatively the validity of the enlire model . Rather, it may be neces- sary to assess individually the correctness ·.of the representation of various mechanisms which together influence the system . This separation virtually necessitates that the overa l l ev~l ~ation will be qualitative . Whereas in assessing model validity the emphasis is placed on t he segments that comprise the model, in assessing model accuracy the emphasis 106. Figure 1 Evaluation of Urban Scale , Photochemical Air Quality Models Fundamental Model Development Scientific Principles Fundamental Correctness Computational Model Development Ability to Replicate Known Conditions Observations Computational Model Implementation Comparison of Value of Information Out put to Cost Costs Model Applicat jon 107 . shifts to the model as a complete unit, with a shift from evaluation of the beginning of the modeling process, scientific principles, to the end , actual output . Accuracy assessment -measures the ability of the model to replicate known conditions . The main limitation to a complete study of accuracy is the availability of data to set the standard against which the model is to be judged. In Section 2 we discuss in detail the quan- titative methods for accuracy assessment. However, even here, it is difficult to arrive at a single quantitative index of model accuarcy. We will always be left to form a qualitative summation of the results of the various methods. As air quality models receive increasing attention in policy decision applitations, the interest in their efficiency also increases. Because of limited resources it is usually necessary to compare the value of a model's information output to the costs associated with running the model . * Effi - ciency may at first appear to be a quantitative measure for model evalu ation; the qualitative aspect enters into the 'value of the information output' part . A model does not necessarily have twice the information content if it has four vertical grid layers instead of two. Any measu re of efficiency ultimately reflects an opinion in defining the important outputs and the constraining costs. · In concluding this discussion of the evaluation framework, we have *We define cost in a broad sense that includes time and resources i.n addition to direct dollar outlays. Thus , t h~obvious cost is in computer expenses, but if three or four months ari required to acquire and assemble the input data base, then this is also a significant cost. Of course , we could assign a do 11 ar value to time and resources and then use a s·ingl e unit of cost. Such an approach becomes artificial, however , when we ar e constrained by time or computerequipment limitations. 108 . two general comments . The first relates to the process of evaluati.on and the second to the results of evaluation. Air quality model evaluation has been presented in this discussion as a serial process, progressing through stages from development to application. As a conceptualization, this allows us to identify and isolate the important aspects of model evaluation. In an ongoing model development effort, the evaluation process will not follow this form explicitly. The various ·stages will often occur in parallel and with ~ considerable feedback . In practice, all of the stages are interdependent. The framework we have presented must be used as a framework and not a fl ows.heet. The result of a model evaluation is not simply a single index of total performance. for this . It is clear that air quality models are too complex In the final analysis, subjective judgements , based on knowledge of the physical problem being studied, are necessary . 109. 2. ACCURACY ASSESSMENT 2.1 Introduction As defined in Section 1, accuracy refers to the ability of a model to replicate known conditions; In some situations it may be possible, due to characteristics of the process or of the model, to specify a single index of accuracy; this is decidely not the case for air quality models . With air quality models there is typically a very large number of predictions and observations. Thus, a given model is generally not uniformly 'good' or 'bad' in all of its predictions. Rather than in- sisting on a one-to-one comparison of predictions and observations, it is important to identify the characteristics of a 'good' model and to assess model accuracy on the basis of these characteristics. In Section 2.2 we will identify a variety of methods which, collectively, form the basis for assessing the accuracy of air quality models. In this work we discuss the direct comparison of model predictions to physical observations. on such comparisons. Accuracy assessment need not be solely based Comparisons can be made relative to input parameters. For example, in a study simulating many days, an index of direct prediction to observation comparison can be analyzed as a function of daily emissions. A model may not necessafily perform as well at one level of emissions as another, and this type of analysis is sensitive to such behavi.or. Mass fluxes also reflect very basic characteristics. of the air qual Hy in an airshed. Observed and predioted mass fluxes can be compared I • by forming the new variable, mass flux, based on the input parameters, emissions and transport variables. 110 . Direct comparisons of model predictions with physical observations are the basic accuracy assessment methods and thus we focus on them. The methods involving the input parameters are important but are not as general nor as objective. Computer codes for the direct comparison can be written independently of det~iled knowledge of the model to be evaluated. Codes involving the input parameters tend to be more model specific because different models, even of the same physical process, will require different details in their respective inputs. Also, the farther removed from the raw data, the more subjective the judgments that may have to be built i.nto the analysis become . . If these judgments involve issues similar to t hose encountered in the model development itself , then the risk exists that the model evaluation is biased by .our judgments . An example is in the determination of the wind Held from the data at a few observation locations. Different models treat this problem differently and so would different evaluation schemes. As mentioned earlier, air quality models are sufficiently complex .._ that accuracy must be assessed by a variety of methods. Before presenting s~ecific methods, we discuss the two levels at which predictions can be directly compared to observations. At the second level the evaluations become ~ore meaningful because wear~ adding to the assessment more implicit information about the nature of the data . Level 1 assessment does not employ any i-j nformation concerning the ' process .being modeled and considers the obse~~ations and predictions to I .j ,' be ;two sets of data which sho~ld inherently .be equivalent. In essence, . at .this level we measure the degree to which the two sets are in fact equjvalent. Level 1 evaluation~ are generally meaningful to perform on 111. models of any process and -they represent the quality of the model's fit to the observations on an average or overall basis. These two character- istics are, respectively, the strength and the weakness on Level 1 methods. Level 2 methods are similar to those of Level 1 in that it is assumed that the two sets of data should inherently be equivalent; how~ver, they do account for the particular process being modeled and the use to which the model results will be put. These considerations lead to methods that consider the aspects of the data which indicate whether or not the model will be useful for its intended application. Thus, in contrast to Level 1 methods, these methods are process specific and only measure the model's performance based on specific aspects of the data. In the broadest sense, air quality models may be used either for (1) evaluating proposed emission control strategies, or (2} increasing our knowledge in the science of atmospheric pollution. current or In this discussion we will be con- cerned with only the first use of air quality models. ft is for this use that comparisons to physical observations are most relevant . Various indices of air quality are used to evaluate the effect of control strategies. Thus, a useful model must produce similar values to the observations for the air quality indices. The Level 2 methods of comparison measure this characteristic of the model. It is extremely important to realize ·that there are three fundamental . i difficulties in comparing air quality observations to predictions of the type we are considering. First, on the sc~le of the model , the observations ( J are spatially poi ntwi se, whereas the predi ct,ions represent volume averages. Second, the observations contain instrument errors . i nput parameters present the final problem. Er rors in the model Even if the model is an idea l 112. fonnulation of the air pollution process, the predictions will be in error if the inputs are in error. In evaluating the results of an accuracy assessment study these factors must be remembered. The presence of these factors is one reason why many diverse methods are needed for a meaningful assessment. Next, we discuss methods for Level 1 and Level 2 assessment. 2.2 Assessment Levels 1 and 2 The list of possible methods for assessing model accuracy is indeed lengthy, the fundamental reason for which is simply the sheer number of data values that must be dealt with. (A relatively modest evaluation effort might involve hundreds of values, and in a comprehensive regional model comparison study the number may reach into the millions.) Thus, the number of potentially meaningful ways to organize and reorganize the information for inspection is definitely substantial. The data could be viewed temporally, spatially, spectrally, or relative to frequency distributions. Superposed on these classes of organization are the many possible forms of averaging. In Table 2 we catalog many of the accuracy assessment methods that have been purposed for comparing air quality model predictions to physica l observations. No claim is made that this catalog is necessarily complete ; l rather, it is representative and includes the' important concepts . From this catalog a general accuracy assessment package , AQMAAP, has been developed. We now discuss the methods tha,t comprise AQMAAP . Included in AQMAAP Quadratic Measure of Mean Error Included in AQMAAP Included in AQMAAP Included in AQMAAP Included in AQMAAP Included in AQMAAP Included in AQMAAP RMS Error Centered at 0 RMS Error_, Cent€Ped ·at Mean ,Errqr - Absolute Error Band Percentage Error Band Residual Histogram Residual vs . Time Plot Residual vs . Observed Magnitude Plot Corrnnents Mean Error Anal1sis of Residuals Level 1 Method Accuracy Assessment Methods Table 2 Liu et al . (1976) Liu et al. (1976) Liu et al. (1976) (1977) Maldonado and Bullin Liu et al. (1976) Johnson et al. (1976) Maldonado and Bullin (1977) Anderson et al. (1977) Maldonado and Bullin (1977) Reference N N (:,,J Included in AQMAAP Residual vs. Location Plot Chi-Square Test Other Anal_,lses Substantial uncertainty in assuming a frequency distribution for the errors. .(Same as comment above.) Numerical Correlogram . . .- Graphic presentation of information found from correlation coefficient and from regression analyses. ,_ Scatter ~lot ~ Included in AQMAAP Regression Analysis ......,,_ Included in AQMAAP Correlation Coefficient Analxsis of Trends Corrments Yields similar information to the residual vs. observed magnitude plot. Method Residual vs. Predicted Magnitude Plot Table 2 (continued) Liu et al. (1976) Anderson et al. (1977) Liu et al. (1976) (1977) Johnson et al . (1976) Maldonado and Bullin Nappo (1974) Liu et al. (1976) Johnson et al. (1976) Reference Liu et al. (1976) N N >I:,. -- Included in AQMAAP Included in AQMAAP Included in AQMAAP Included in AQMAAP Included in AQMAAP Magnitude and Time of Concentration Peaks Location of Concentration Hotspots Exposure Geometric Parameters of Highest Values Frequency Distribution Analysis -~f Indices of Air Quality ~ Requires data over periods long enough to have many cycles. A model accuracy study will typically have data for only a few days. Spectral Analysis ....,. Yields information similar to methods for analysis of residuals. Fractional Mean Deviation from Perfect Correlation Level 2 Comments Yields similar information . to methods for analysis of residuals. Method Mean and Variance of Ratio of Predictions to Observations Table 2 (continued) Anderson et al. (1977) Nappo (1974) Reference ~ ,_, ,_, 116. The basic nomenclature for the concentration values is Observed concentration of specie i, at location x., J on day dk and at time t 1 . ei(xj,dk,t 1): Predicted concentration. Mi: Number of species . Mj: Number of locations. Mk: Number of days. M1 : Number of hours of interest during each day. Many of the methods will represent averaged properties of the data. Because there are three independent indices (location, day, hour) there are many ways to treat the data in determining averaged properties. such ways appear particularly useful for air quality data. to consider the data in Five It is important this many ways to first have a general measure of the model performance and then to look at more specific measures. In this manner, systematic or important isolated deficiencies in the model may be uncovered. Most of the Level 1 methods in AQMAAP will treat the data in each of the five ways. To avoid needless repetition we use a notation that allows each of the five to be denoted by one set of equations. lists Table 3 the five averaging processes along with the interpretations of the shorthand notation. The meaning of the notation will become clear once we have discussed the measures of accuracy assessment . • The Level 1 methods are classified int_o two types: residual and analysis of trends. analysis of In the fitrs,t, we ex·amine the differences between the observations and predictions. Ideally, the residuals would all be identically zero. In the second, we conceptually order the observations from smallest to largest, then examine the degree to which the (A) Measure Notation: Averaging operator. Function of the sets of indices m and n. Set of indices to be averaged over. Set of indices si is a function of . ~L m=l g m n M. Function of the averaged function g. J g(m,n~ f m=l t ~ I~ Measures for species i and set of indices n. = f si(n) si(n) Averaging Shorthand Notation Table 3 h-1 h-1 "-J Location Measure over all locations. 1 ~I: J j=l J 1 Day Measure' over each hour. Ml M L l l=l Mk 1 z:~1 Hour Measure over each day. L L k k=l l l =l 1 1 M LM M L k k=l M. H z:i ML 1 M L m=l 1 fM. 1 Mk 1 M1 Mj Mk M, j=l k=l l=l Mk M1 z:~1 s'.1 Day, hour Measure over all times. Symbol of Measure s:L 1 Indices Averaged (Set 'm'} Location, day, hour Averaging Notation Measure over all time and locations. (B) Table 3 (continued) k = l, ... ,Mk;l = l, ... ,M 1 j = 1, ... ,Mj; l = 1, ... ,M 1 j = 1, ... ,Mj;k = l, ... ,Mk j = l, ... ,Mj Set 1 n 1 Co h-..1 h-..1 119. corresponding predictions follow the identical ordering. (A) Analysis of Residuals Denote the residual as wi where Mean Error Mean error measures the average bias in the predictions and indicates whether the model predominantly over or under predicts. form of the mean error is The general Root Mean Square Error Centered at the Mean , This RMS function measures the average spread of the residuals. It is insensitive to the bias. o.(n) l The general form is defined = *As an example of the shorthand notation, g is the residual w. and there 1 is no function f. The mean error over all times is Mk M1 ii'.(x.) 1 = J J-k k=l ~ J- L wi(xj,d'k,t l l=l · 1) and the mean error over all locations is M. ~~(dk,tl) = J. ~ wi(xj,dk,tl) • J j=l In this equation, f is the square root a~d] is the squared difference between the residual and the corresoonding mean error. Thus, M· · -rq · M • 1 ** I a;TL = ~ ~ 2- M1 ) k M1 l (w; (xj ,dk, t 1 ) - oiTL) 2]'2 • J j=l k ti. l l =l [ 1 L 1 20. Error Bands In contrast to the above averaged measure of the spread, we can evaluate the model relative to a prescribed degree of spread by finding the percentage of the predictions which fall within a given band around the observations. ± yi This band may be based on an absolute concentration, ppm, or a percentage, ± zi % of the observed value. The absolute error band measure is ai(n) = the percentage of residuals from set 'm' The percentage error band measure is s (n) = the per1 centage of residuals from set 'm' within ± zi% of ci(m,n). within ± yi ppm. Residual Plots Residual plots provide a summary picture of the nature of the bias and the spread for the model . value for air quality models . AQMAAP: Five types of plots are of specific Two residual histograms are included in a histogram based on all of the residuals for a given specie, and individual histograms for each location. By plotting the residuals in other ways we can determine if there are certain condi'tions under which the model is good, while others under which it is poor . Thus , we also plot the residuals versus time of day, versus location and versus observation magnitude . (B) Analysis of Trends Methods for analyzing resi.duals pf f.m,arily indicate the absolute I fit of the predictions to the observations : • ' Predictions based on funda- ,. mentally wrong relationships can have measures of the residuals that appear similar to those for soundly based predictions . By analyzing t~e 121. trends of the predictions relative to the observations we gain further information indicating whether or not the predictions at least obey the same relationships as the observations. Correlation Coefficient The correlation coefficient measures the degree to which the magnitude of the predictions increases linearly with the magnitude of the observations . The coefficient, however, is insensitive to the extent of the increase. If the predictions increase linearly at only 1/10 the rate of the observations the correlation coefficient will still be one. The realtionships defining the correlation coefficient are 1 M pi(n) = M M ~lL Mm=l v.(m,n) l . L [v;(m,n)ni (m,n)] m=l L (ni(m,n)} {vi (m,n)) 2 ~ m=l = ei(m , n) - i\(n) ei(n) = ~ 2 ! '2 ] M I: e,(m,n) m=l ni(m,n) = ci(m,n) - c:; (n) c:; (n) = 1 M M I:ci(m,n) . m=l Linear Least-Squares Curve Fit l In addition to the correlation coefficient, we can also measure the average increase in the observations a~ . the observations increase . I 1 ,' The slope parameter of the linear least -squares curve fit is th i s measure . The intercept parameter measures the bias if the sl ope parameter is very nearly one. Otherwise , the ~lope so strongly influences the intercept 1 22. that it has no particular independent interpretation. We seek to minimize when ·~ is defined = e.(n)c.(m,n) + ¢.(n). l l l Thus, the slope and intercept are 1 M ML[v;(m,n)ni(m,n)J m=l = 1 M M I:(v;(m,n) )2 m=l and ¢i(n) = ~i(n) ei(n)§i(n) , respectively. Aside from the plots of residuals tihe above methods can all be applied in each of the five manners of averaging in Table 3. Of course, we can now further average the quantitative measures themselves. Par- ticularly popular are spatial averages of the temporal averaged measures and the converse . The value of such repeated averaging is, howeve r, 'l dubious because physical interpretattons of ·the measures become tenuous and too much averaging tends to obscure the v-ery deficiencies one i s 1 ' seeking to identify. ( The Level 2 methods assess a model relative to its :potential use for evaluating control strategies. These methods measure the ability 123. of the model to produce the same values as the observations for air quality indices. Analysis of Indices -of Air Quality The most common indices of air quality relate to maximum values. We consider spatial and temporal maximum values. Concentration Hotspots On an hourly basis we determine the locations of the observed and predicted hotspots. Then we report yi' the percentage of times the model accurately predicted the location of the hotspot of specie i. Concentration Peaks At each location we determine both the magnitude and the time concentration peaks. The predictions may be basically accurate near the peak but still fail to reproduce exactly the peak value itself. Instead of basing our judgment on only one value we also characterize the group of highest values for a specie. Because air quality data may span more than one order of magnitude, geometric parameters are appropriate. We compare the geometric mean and standard geometric deviation of the upper ten percent of the data. The peak measure~ are 1 magnitude error, £/xj ,dk) " (mtx ~/xj ,d~ ~i:])] - mtx(C\ (xj'dk, t 1 1 • 1 time error, T.(x.,dk) l J = ([ttme of c max] - &ime of C maxJ} l]} 124. geometric mean, MIO A. (x.) 1 J = exp ft/ L ln c.(x.,m)'l 1 J 'J l_l 10 m=l MlO ~- (x.) 1 J = exp [Ml L ln e.(x.,m0'J l_i 10 m=l 1 J and standard geometric deviation, 1/J.(x.) l J = exp{~l 10 A 1/J. ( X . ) = ex{~~o MIO = (MkM, )/10 l where J Ji.:} MlO ~ ( lnci(xj,m) - ln\(xj) ) 2 2 • ~k} 2 MlO ( . lne.(x . ,m) - ln~-(x.)) 2 m= 1 1 J l J L Frequency distribution Another standard manner of indexing air quality data aside from maximum values is relative to their frequency distributions . It is common , for example, to seek the concentration that occurs at the 99th, 90th or 50th percentile and, conversely, to seek t~e percentage of time that certain air quality levels are exceeded. ~ll of this information is compactly summarized on a frequency distri but ion plot. These plots are drawn for each location. Exposure The final index of air quality that we will compare is a simple 125. measure of exposure. We compare the ppm-hour exposure of a receptor at each of the observation locations . o.l (X.) J = 1 Mk 8.l (X.J ) = 1 Mk Mk ~ Ml Ml L Ci ( j 'd k'tl ) • ( tl - tl -1) k= 1 l 1 =1 Ml Mk 1 X L <\(xj,dk,tl)·(tl - tl-1). L ~ k=l l ==I 126 . 3. SUMMARY The final result of this work is the development of the FORTRAN code AQMAAP. AQMAAP represents a general and extensive package for assessing the accuracy of air quality models. The discussion of the total model evaluation framework places AQMAAP in the proper perspective for interpretation of its output. The methods presented here are both meaningful and practical for studies ranging from small one model validations with a few days data to large multiple model comparisons with possibly a few months data. package is also quite general and objective. The Essentially, it requires only the observed and predicted concentration values for inputs . No information on the model or details of the field study need to be provided. For these reasons, any user can assess any air quality model without extensive alterations and with a minimum of subjective judgments. 127 . REFERENCES Anderson G. E., Hayes S. R., Hillyer M. J., Kill us J. P. and Mundkur P. V. (1977) Air Quality in the : Denver Metropolitan Region 1974-2000, Environmental Protection Agency Report EPA 908/1-77-002. Brier G~ W. (1975) Statistical ·uestions Relatin to the Val i dation of Air , Qualito/ Simulation odels, Environmental Protection Agency Report EPA 650/4- 5-010. Johnson W. B., Sklarew R. C. and ·Turner D. B. (1976) Urban air quality s.imulation modeling. Air P01lution, Volume I (Edited by Stern A.C.) . Academic Press, New York, 503~562. Liu M. K. , Whitney D. C., Seinfeld J. H. and Roth P. M. (1976} Continued Research in Mesoscale Air Pollution Simulation Modeling : An~lysis of Prior Evaluation Studies, Environmental Protection Agency Report EPA 600/4-76-016a. Maldonado C. and Bullin J. A. (1977) Modeling carbon monoxide dispersion from roadways . Environ . Sci : :rechnol. 11_, 1071-1076. Nappo C. J., Jr . (1974) A method for evaluating the accuracy of air pollution prediction models. · ·Preprints of the Symposium on Atmospheric Diffusion and Pollution, Santa Barbara, CA, American Meteorological Society, 325-329. ' 128. APPENDIX. A.1 AQMAAP, Air Quality Model Accuracy Assessment Package Description AQMAAP is a FORTRAN coded package of methods for the assessment of the accuracy of air quality models . The scope and attributes of the principles behind AQMAAP have been discussed in the preceding sections of this chapter. This appendix describes the FORTRAN package itself. AQMAAP consists of 41 subroutines comprising the individual methods, 15 supporting subroutines, and 7 main or driver routines. In Table A. 1, the assessment methods and the corresponding subroutines are listed. The package is running on the IBM 370 at the California Institute of Technology Computing Center . Aside from bookkeeping infor- mation (i.e., identification of the data sets and specification of the methods to be exercised), the only inputs are the observation and prediction data sets. A few statements (mainly a single group of DIMENSION statements) in the code are dependent on the size of the data sets . In the following section the necessary information for operating AQMAAP at any installation for any size data set is discussed. 129 . Table A.1 AQMAAP Asse.ssment Method Subroutines Assessment Method LEVEL 1 Mean Error Over all times and locations Over all times Over each day Over each hour Over all locations Root Mean Square Error Centered at the Mean Over all times and locations Over all times Over each day Over each hour Over all locations Subroutine Name Elllll EllOll El2011 El3011 ElOlll El1112 E11012 E12012 E13012 El0112 Absolute Error Bands Over all times and locations Over all times Over each day Over each hour Over all locations E11121 Ell021 E12021 El3021 E10121 Percentage Error Bands Over all times and locations Over all times Over each day Over each hour Over all locations E11122 Ell022 E12022 E13022 E10122 Residual Histogram Over all times and locations Over all times Elll31 El 1031 130. Table A.1 (continued) Assessment Method Subroutine Name Residual Plots Residuals versus observation magnitude Residuals versus location Residuals versus time El 1133 Elll34 [11135 Correlation Coefficient Over all times and locations Over a11 times Over each day Over each hour Over all locations E11141 [11041 E12041 El3041 El0141 Linear Least Squares Over all times and locations Over a11 times Over each day Over each hour Over all locations El 1142 Ell042 El2042 El3042 El0142 LEVEL 2 Magnitude of Daily Concentration Peaks Time of Daily Concentration Peaks Location of Concentration Hotsoots Peak Geometric Parameters Daily Exnosure Frequency Distribution Plot. Over all times E22051 [22052 [20156 E21062 E22071 [21061 131. A. 2 Operational Aspects Tables A.2 - A.5 detail the operation of AQMAAP. Tables A.6 and A.7 provide additional information that may be pertinent for operating AQMAAP at another installation. The operation of AQMAAP requires that a few parameters be set by the user. Most of these relate to the size of the data set. The remaining parameters give the user control over the error band criteria and some flexibility in the presentation of the plots. In Table A.2 the program parameters are listed. The three inputs to AQMAAP are, first, the identification of the data sets, then, the data sets themselves and finally, the specification of the methods to be exercised. The proper input sequence and FORTRAN fonnats are listed in Table A.3. Five averaging processes are used in AQMAAP. As discussed in Table 2.2, a different set of indices is averaged for each averaging process. The corresponding measure of accuracy is then computed as a function of the remaining indices {including the specie index}. These indices are specified in Input section {3}. Any partic- ular method need only be exercised over a selected portion of the data set. For example, the mean error evaluation over all times might only be computed at a few locations instead of at all locations . Table A.4 elaborates on these input codes. A sampl~ input deck is shown in Table A.5. Table A.6 indicates which of the 15 ~upporting subroutines are used I in each of the main routines. • ' Eleven addittonal routines available at the California Institute of Technology Computing Center {CITCC) are used in AQMAAP. These routines are referenced and cataloged in Tabl e A. 7. (3) (2) (1) ..... ...- DIMENSION {~:EQ\ (Mc), X(Mc), Y(Mc) RES) where M = M.MkMl and NR is t he vari able used in Elll33, FREQ in E21061 and RES in the remain i ~g fotlr subroutines. . DIMENSION WS2R(M1 ,Mk,Mj,Mi) In the subroutines E11131, E11031, E11133, E11134, E11135 and E21061, specify the dimensions in the statements where Mb is the greater of Mi and 3. COMMON--, /DATA/ COBS(Mi ,Mj ,Mk ,M1 ), CPRO(Mi ,Mj ,Mk ,M1 ) COMMON /WRKSPC/ WS 1(Mb ,Mj ,Mk, Ml ) , WS2 (Mb ,Mj , Mk, M1 ) COMMON /CODE/ NNSPC(Mi), NSSPC(Mi), NSPC(7 ,Mi), NNLOC(2 ,Mj), NLOC(7 ,Mi) , NALOC(2 ,Mj), NDAY(7 ,Mk) ~ In the ma i n routines, the method subroutines and support subroutines INCOD3, INDATl, OTDAT3, RESID~nd OTRES3 specify the dimensions in the statements where Mi i s the number of species (in the data set), Mj is the number of locations, Mk is the number of days , M1 i.s the number of hours and Ma = Mi M/1kMl . INTEGER I PARM/ Mi /, JPARM/ Mj /, KPARM/ Mk /, LPARM/ Ml/ REAL COBS(Mi ,Mj ,Mk,Ml )/t\ *-1.0/, CPRD(M; ,Mj ,Mk,Ml )/Ma *-1.0/ In the BLOCK DATA section of the main routines, specify the values and dimension in the statements Data set dimensions Program Parameters AQMAAP Table A.2 N C"1 I:\:, In E11134 (residual versus location plot) and Elll35 (residual versus time plot), the length of the abscissas equal the number of locations and days, respectively. These subroutines are currently written to plot on a single sheet of paper. If there are more than 55 locations or 3 days, it will be necessary to make alterations to either the plotting format or the available plotting area. Such changes will require additional familiarity with both the plotter being used and with the coding of these two subroutines. ...:._,._ .,... . ,.._~~ where yi and z1 are, for each specie, the user selected sizes of the absolute and percentage error bands, tespectively. EBP(i) = zi' i=l, ... , i\ EBA ( i ) = y i' i =1 , ... , M; In EB, Exxx21, Exxx22 and SETEB, specify the dimensions in the statement COMMON /ERRBND/ EBA(Mi), EBP(Mi). In SETEB, specify the values in the statements In Elll33 , Elll34 and E11135, specify the values in the statement s MSB(m) = b(m), m = 2,3 ,4 MSE(m) = e(m), m = 1,2,3 with values between 1 and 99. Four symbols are available to these plotting routines to indica t e the frequency of a given residual value. The mth symbol is drawn if b(m) ~ (the f requency of the residual) ~ e(m). Plotting, sy~bol specification (2) (1) Error band size (4) Table A.2 (continued) ('.>J . hl l'.,;i where v· is the (plotting) position along the abscissa in E11134 of the j th (observation) locatio~. This allows the user the option of grouping the locations in some meaningful manner (e.g. along a typical wind trajectory or ordered in distance from a central point). MLOC(j) = vj In POSIT, specify the value in the statement Abscissa QOSition Table A.2 (continued) N <:,J >t:,,. 1 513, 212, IX, 6El0.3 5I3, 212, IX, 6El0.3 lOX, 7A4, 14X, I4 Day information Oay' ',-"-i ndex- code number (2) Data sets Observations Specie, Location , Day, Hour start , Hour finish, Units , Method, Concentrations (ppm) End of data flag; •~-1• Predictions Specie , Location, Day, Hour start, Hour finish, Units, Method , Concentrations {ppm) End of data flag; ~-1' 12, IX, I5, 2X, 7A4, 14X, I4, 4X, 2A4 I2, 2X, 1A4, 2X 7A4, 14X, I4 Format Location information Agency code numbers (County and Site), Name, Index code number, Abbreviation In,eut (1) Identification of the data sets Specie information Agency code number , Symbol, Name, Index code number AQMAAP Input Sequence and Format Table A.3 1 Dependent on amount of data used. 1 Dependent on amount of data used. KPARM JPARM IPARM Number of Cards <:Ji ('.>J ....... lnQut (3) Specification of methods to be exercised Method subroutine code For this method: Number of species for which method will be exercised Specie For this specie: Number of locations for which method will be exercised Location For~ h_is sp~cie/location: , Nu~ber of days for which method will be exercised Day For this location/day: Number of hours for which method will be exercised Hour Table A.3 (continued) 40X,I5 45X,15 30X,15 35X,I5 20X,I5 25X, 15 lOX ,15 15X,I5 IlO Fonnat ~ Dependent on output requested 1 Number of Cards >--, c,;i CT;, To denote a missing concentration value enter 1 -0.lOOE+Ol 1 • a str ing of missing concentrations. (6) (9) To skip a method, specify 'zero' for the number of species. For Input section (3), the required inputs vary with the specific method. (8) See Table A.4. Compl~te Input section (3) for each method subroutine in the program in the order given i n Ta6l e A.4 . (7) ,- The assumed units ~or all concentrations are 'ppm'. (5) - ... Use the appropriate index code numbers for designating specie, location and day in Input sections (2) and (3). (4) .. . In Input section (1),agency code _number, units and method are optional inputs. (3) ~ In Input section (1) the index code numbers are chosen by the user. They must be unique with 1 ~ i ~ IPARM, 1-~ j ~ JPARM and 1 ~ k ~ KPARM. (2) Delete the entire card for A complete input deck is required for each of the 7 main routines. (1) Notes: Table A.3 (continued) ,__, 'J <J-J 138. Table A.4 AQMAAP (A) Sequence of method subroutines to be followed for Input section (3) Main Routine MN - Basic Measures of Bias and Spread Subroutine Sequence within Main Routine Elllll E11112 EllOll E11012 E12011 E12012 E13011 E13012 El0lll E10112 EB - Error Bands E11121 E11122 E11021 El1022 E12021 E12022 E13021 E13022 [10121 E10122 HT - Residual Histograms E11131 Ell031 RF - Residual Plots Ell133 El1134 El 1135 CR - Analysis of Trends E11141 E11142 E11041 E11042 E12041 E12042 E13041 E13042 E10141 E10142 139. Table A. 4 (continued) Subroutine Sequence within Main Routine E22051 E22052 E20156 E22071 E21061 Main Routine IN - Indices of Air Quality FD - Frequency Distribution Plot (8) Sequence of the data set indices to be followed for Input section (3) Type of Evaluation Over All Times and Locations Over All Times Over Each Day Over Each Hour Over all Locations Method Subroutine Code xllxx xlOxx x20xx x30xx xOlxx 1 1 1 1 1 2 2 2 Specie X Cl) E21061 Location "'Cl C: ...... Day Hour 3 2 3 3 140 . Table A. 5 AQMAAP Sample Input Deck for Program 'MN' (1) - 19 PTA POLL~TANT A 0~01} 46 91 POLLUTANT D PTO 66001 LOCATION 22 0004 00_01 91 75999 LN 22} LOCATION xx 0003} LN XX 11 MONTH N 0000 0001 22 MONTH N 0000 0002 1 05 1 1 O.OOOE+OO 0. 040E+Ol 0 .080E+Ol O. lfOE+Ol 0.160E+Ol 0.200E+Ol} or ( 2) - 4 -1 3 3 2 18 23 1 2 0.200E+02 0 . 100E+02 !.OOOE+02 0 . 900E+02 0 .800E+02 0 . 900E+02 ~ 11 1 0 0. 760E+00-0.100E+Ol 0. 520E+OO 0. 6?0E+00-0.100E+Ol 0.420E+OO} -1 (3) - Species, Pollutant A is index #1. Locations, Location XX is index #3. Days, The 22th is index 12. Observations ( as available for 4 species, 3 locations, 2 days and 24 hours ) . END OF DATA Pred i ctions (Note : No data for hours 0700 and 1000). END OF DATA 11111 3 Exercise method El 1111 for the 3 species , 1,2 and 4. 1 2 4 11112 0 Skip methods Ellll2, EllOll and El! Ol2 . 11011 0 11012 0 12011 2 4 2 Exercise method El2011 for: spec ie 4, location 1, day 2, specie 4, location 1, day 1, specie 4, location 3, day 1 and specie 3, location 2, day 2. 2 2 1 3 3 2 2 } 12012 0 13011 0 13012 0 Skip methods E!201 2, E!3011 and El301 2. 10111 2 2 6 06 07 08 09 lj) 11 Exercise method E!Ol 11 for : specie 2, day 1 and 0600 hours -+ 1100 hours and specie 2 , day 2 and 1200 hours -+ 1700 hours 2 6 10112 \ 12 13 14 15 16 17 Skip method El0112 . 141 . Table A.6 AQMAAP SUPPORT-MAIN REFERENCE TABLE MAIN ROUTINE INCOD3 . INDATl OTDAT3 RESID OTRES3 l.J..J z ...... SETEB I:::> POSIT 0 0:: co AXES :::> TICKS I0:: 0 a. SCLSl a. :::> SCLS2 SCLS4 SCLS5 SCLS6 SCLS7 MN EB HT RF CR IN FD X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X V') V') X X X X X 142. Table A.7 AQMAAP {A) AQMAAP-CITCC REFERENCE TABLE El1131 E11031 Ell 134 El1135 E11133 E21061 AXES TICKS SCLSl SCLS2 SCLS4 SCLS5 SCLS6 SCLS7 Mfl LLJ EB z HT ..... RF I-:::, 0 0::: Cl. c:i:: c:i:: ~ 0c:i:: CR IN FD E21062 Exxxll Exxx12 Exxx21 Exxx22 Exxx41 Exxx42 E22051 E22052 E20156 E22071 INC0D3 INDATl OTDAT3 RESID OTRES3 SETEB POSIT X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 1 43 . Table A.7 (continued) (B) Catalog of CITCC Routines Routine (Name/Number) RTIME* C 868-249-370 Compute remaining execution time. SORTDE Cl268-279-370 Sort an integer array in decreasing order. SORTI2 Cl268-279-370 Sort and index a real array in increasing order. SORTIR Cl268-279-370 Sort a real array in increasing order. YNTERP C 766-176-370-10 Evaluation of functions available in tabular fonn at arbitrary values of the independent variable. MAXMIN C 267-214-370 Detennine the maximum and minimum of an array. OUTCOR C 169-288-370 Convert data with a FORMAT and place the EBCDIC formated line image in core storage. SYSSYM C 367-218-370 Plot alphanumeric symbols. SYSENO C 467-221-370 Terminate a plot. SYSPLT C 467-222-370 Move plotter pen. XYPLT C 774-411-370 Plot point datq in a 2-0 field. Function *An optional routine, can be removed -without altering basic program function. ,J 144. Chapter 3 On Frequency Distributions of Air Pollutant Concentrations 145. Atmosf'Mr1r En1•iro111nnit Vol. 10, pp 941 - 9.SO. Pcrpmon PrCM 19 76. Pnnttd in Grca1 Brnain. ON FREQUENCY DISTRIBUTIONS OF AIR POLLUTANT CONCENTRATIONS KENNETH E. BENC'ALA and JOHN H. SEINFELo• Department of Chemical Engineering, California Institute of Technology. Pa.<adena. CA 91125, U.S.A. (First rec-ei,,ed 21 Januar)' 1976 and in final form 2 Jun e 1976) Abstract-Observed frequency distributions of air pollutant concentration levels arc critically analyzed with respect to their statistical description. It is demonstrated that several common distributions can be used to fit observed data. one of which is the popular log-normal distribution. The observation that concentration distributions for all averaging times are approximately log~normal can be explained if the short averaging time data are themselves assumed to be log•normally distributed. The near log-normality of pollutant concentration frequency distributions can be ex.plained on the basis of the near log-normality of wind speed distributions. although this explanation does not establish that wind speed distributions are solely responsible for observed concentration distributions. It is concluded tha1 pollutant concentration frequency distributions are the result of complex phenomena and cannot be predicted exactly, but that the approximate log-normal character of the distributions is useful from a practical point of view and can be understood qualitati vely on the basis of the relation between wind speed and concentration. INTRODUcnON Most current United States federal air quality standards are stated in tenns of the yearly frequency of violation of a specified concentration level for a given averaging time. For example, the 1-h average carbon monoxide concentration may only exceed 35 ppm once during the year. This is equivalent to specifying that the 1-h average CO concentration may exceed 35 ppm only 0.01 I% of the time. Therefore, if the frequency distribution for hourly average CO levels in a particular urban area could be predicted, then it could be ascertained what degree of emission control would be required to enable meeting of the air quality standard. There has been interest for some time in the frequency distributions of air pollutant concentrations. Larsen (1971) carried out a comprehensive analysis of the data collected in the Continuous Air Monitoring Program (CAMP) for the years I 962-- 1968. The • To wh om correspondence should be addressed. CAMP data contain measurements of seven pollutants (CO, NO, NO,, NO,. oxidant, SO,, and hydrocarbons! in eight cities. Readings were recorded every five minutes and then averaged over time periods ranging from ten minutes to one year. Based on his analysis of the data Larsen concluded that regardless of pollutant, city, or averaging time urban concentration frequency distributions could be universally represented as the log-nonnal. (The log-nonnal is a two parameter distribution and is represented by a straight line on log-probability paper.) For example, Fig. I shows frequency distributions for one hour averaged CO concentration in various cities. Figure 2 shows frequency distributions for CO concentrations in Chicago for various averaging times. Quantitative explanations of why the distributions tend to be log-normal are incomplete, and there is currently no way of predicting how the distribu1ions will shift if emission levels are changed. The objectives of this work arc as follows . First, we wish to consider the question of wh_ether or not 100.-----,--~-~~~~~---,----,---~~~ o Ctucooo ~ ~ 0 - ~ a: 10 i;;; 5 m,n ov•rooe o J 00r ove,oo• o I t-ovr ove,09e • I month overoge • 8 1\ou• ove r09e o Oenve• • Son F,onc,sco Log-no,mol d1s1r ,buhons using pcrome1e•~ given br Lo•sen U971J. w u z 8 do':-c1---::-0"=10,-~---,1"'=0----'--~~50~~~90~-~99-99~ 90 - 99 _J99 FREQUEN CY, percenl Fig. I. I h average CO concentration disuibutions: CAMP data (1962- 1968). 'o.01 0.10 10 50 90 99 99.90 9999 FREQUENCY, percenl Fig. 2. CO ' concentration distributions in Chicago for varr u? averaging times ; CAMP data (1962- 1968). 146. among common distributions the log-normal does, in fact, provide the best fit for pollutant concentration data. Specifically, we desire to see if other available statistical distributions provide a better fit to a selected sample of data. Second, we wish to analyze the effect of averaging time on the frequcnc-y distribution of air pollutant concentrations. We seek specifically to understand the observation that frequency distributions tend to be log-normal regardless of averaging time. Finally, we desire to analyze the possible physical reasons for the near log-normality of ai r pollutant data. Our overall aim is to attempt to shed some additional light on the issue of air pollutant concentration frequency distributions. From the outset it must be recognized that these distributions are the result of many complex phenomena. and that we cannot expect to be able to predict them exactly in a given situation. Nevertheless, if the extent of validity of the distributions can be understood at least in a semi-quantitative manner, the use of the distributions can be facilitated. REPRESENTATION OF CONCENTRATION DATA Although the log-normal distribution has generally been used to represent air pollution concentration frequency distributions, there are other common distributions which resemble the log-normal and are candidates for representing the data (Lynn. 1974; Mage and Ott, 1975 ; Pollack. I975). It is of interest to examine how well different distributions actually fit air quality data. Because chemical reaction behavior may affect the form of concentration distributions. we confine our attention to CO. Our purpose in this section is to examine the ability of a variety of sta tistical distributions to fit selected air quality data. The statistics of a set of air quality data may be analyzed in terms of either its probability density function (pdf) or its distribution function. The distribution function is an integral function of the pdf. thus for the purpose of evaluating the fit of mathematical forms, comparison of the data to the pdf provides O lhovr CO (1010 lo , Ph ilodelph,o {1 962-1968) o 2-pcromtter roq ·normot dl$trib.1hoo uSf'IQ poromele,s given by lorseo 1197 1) b Leosl-~res bell Ill 3-p()l'Ol'Tl@le, IOQ· normol diSlrtbuhon c Leoi.t·s,qvores ~ , f ,1 QOl'T'fflO drstr1but1on ,o• ,_....,___~_....,__....,__~~....,__~~---'._.__,.__~..., 0.01 0.10 ! 10 50 90 FREQUENCY, percent 99.90 99.99 Fig. 3. CO concentration distributions in Philadelphia with Hoes representing the two-parameter log-normal. the three-parameter log-norm al and the gamma distributions. O lhOUf co dOIO lo, SI Loo1s 11962-1968) o 2-poromeler lOQ·rorrnol drs trobuhoo vSU'll;J pcrometers o,ven by Lo"'en ( 1970 b least· squ:,res bes! I ,I 3 · perometer IOq· normo• d,51r,buhon C Leo1l · SQUOffl\ bell frl pnmo dstrlbution 10 ~'= _0-,-='o.-=,o-~,--,Lo....,__~~~o....,__~-90~~ -,L,-,-,-'-.90 --'99.99 FREOUE NCY, pe,tent Fig. 4. CO concentration dis tributions in St . Louis with lines representing the two-parameter log-normal. the threeparameter log-no rmal and the gamma distributions. a much more stringent test than comparison to the distribution function . Unfortunately, a plot of the pdf is more difficult to interpret than is a plot of the more commonly used distribution function: therefore, we consider the distribution function. In attempting to determine which parameter values in a particular mathematical form provide the best fit of that form to the data one has considerable latitude. We selected an unweighted least-squares criterion for determining the parameters in the distributions which best fit the data. (However, an implicit weighting arises in that relatively more data points were available at high concentration values than at low.) Table I presents four potentially applicable pdfs. Figures 3 and 4 show Philadelphia and St. Louis CO data plotted against three of these distributions, the two-parameter and three-parameter log-normal and the gamma. The parameters used for the two-parameter log-normal plots are those given by Larsen, while for the other two the least-squares best fit parameters are used. On log-probability coordinates the Weibull and gamma distributions are concave while the three-parameter log-normal can be either concave or convex depending on the sign of the third parameter, o. Of the eight data sets investigated, six are to some degree concave while two (Los Angeles and Washington) are slightly convex. Table 2 suggests possible physical interpretations of these deviations from linearity. Table 3 presents a comparison of the sum of squares error in fitting the distribution functions to the eight data sets. In the least-squares sense the three-parameter log-normal is superior to the twoparameter distribufons. The added flexibility afforded by a third parameter accounts for this. In comparing Figs. , and 4 with ihe quantitative measures of goodness of fit given in Table 3. we see that in most cases the two-paramter l; g-pormal distribution provides a useful, if not excellent, approximation to the data but that in that sorhe 'cases it is possible to find other two-parameter dislributions which provide better fits. 147. Ta hie I. Common prohability dcnsily functions Name 2-paramcter x -.. In (,/ r::2n x In /i) _ 1 exp { - (In - -•, . >I' . } 2 ln ·// log-normal J ~ + o ! 3-parameter log-normal -• [ ,_; 2n(I _ { (ln(1+M-In,J ' } b)lnP] 1 exp - -·-····-··, ··-· d1 2 ln ·p r.= . lv ,rn (x + b) In Pl 1 exp { - x=O 2 (In Ix + b) - In >1 } --· ·-i"J;:Zf-·· x> O X > - b Weibull Gamma Table 2. Implication of the shape of a distribution function represented on log-probability coordinates relative to a straighl line (the log-normal} Concave shape At high concentrations Convex shape Fewer days of high concentration. More day s of higher concentration. An upper limit in concentrati on suggested. At low concentrations. More days of low concentration. Fewer day s of low concentration. A lower limit (back~round) in concentration sugtz:ested. Given the complex dynamic nature of urban air pollution, no one distribu1ion will always be the best or even an adequate representation of the data. however, the two-parameter log-normal distribution is clearly a useful distribution both for describing data and for establishing an understanding of air pollutant statistics. EFFECT OF AVERAGING TIME ON FREQUENCY DISTRIBUTIONS Ambient data are generally reported in terms of time averaged values. The time averaged value of concentration C centered at time r and over an averaging period T can be represented as If"'' C(~)d~. <'.',{r) = - Typical values of T in air pollution applications range from 5 min to I y. It has been observed that air pollutant concentration distributions approximate the log-normal regardless of averaging time. and that the median concentration is proportional to the averaging time raised to a power (Larsen. 1971 ). Based on this empirical observation. the standard geometric deviation a,, for one averaging time T, can be related to the standard geometric deviation a, . for another averaging time T. by (2) where V = (l_n(T/ T,l)' '' . (I) ln(T/ T0 ) T , - T/ 2 Table 3. Sum of squares error in fitting ·the distribution s in Table I to I-h average CO CAMP data. 1%2- 1968* Two-parameter log-normal Larsen (1971) City Los Angeles Philadelphia Denver San Francisco Cincinnati St. Louis Washington Chicago values Best fit 0.35 0.15 0.87 0.76 0.64 1.25 0.31 7.24 0.12 • Error based on reduced variate. O.Q7 0.20 0.56 032 0.44 0.08 1.1 7 Three-parameter log-no, mal Weibull Gamma 0.03 0.0 1 1.08 0.55 O.Q3 OJ6 0.14 0.3t 030 1.14 0.13 0.78 0.08 0.48 0.22 0.20 0.17 0.8 3 0.04 0.57 0.20 'o.04 0.05 0.04 (3) 14 8 . where T is the total period over which data are available (usually one year). 50 It is of interest to examine if this empirical observa,. 09 . ,, tion can be explained strict ly on the basis of the - Ft' properties of log-n ormali1y distributed random vari10 ables. Thus, we ask--if the raw data averaged over a period of. say, 5 min are assumed to be log-normally distribu ted, wi ll the data averaged over longer periods continue to be log-normally distri buted. Let X(t) represen t the average value of the concentration over the time period from r - T to r. Thus. ci_o;;;1---,,o." -10-7 --,,;;_o~~-c-,50~~.....,,90,,....-ct99~99""" _90~99=-!_99 X(t) is considered as the ··raw" data. with an inherent FHEOUENCY Ot STRIB UTION FUNCTION , percent averaging time T. attributable to instrument function. The record of raw data then can be represented as Fig. 5. Comparison of F1 to F~ for a typically low va lue.' of a,_,. the sequence. X(I,). X(t 2 )• ...• X (1.). where X(li) is the value in the interval [r, .r, + , ]. X(/ 2 ) is the value Using (5- 7). we obtain Fz as in the interval [r, + ,. r, + 2T] = [1,.r , + ,]. etc. Now. the average concentration over the double interval r - 2, to r is written as Z(r) = Fz<: ; r) 2. ,_, 2 ,r In x·' ½[X(I - T) + X(I )]. The series 2(1 2 1.... , Z(t,) then represents the sequence of concentrations averaged -(lnx - Inµ ,., )') x ex p ( 21n 'a,, - • over periods of length 2,. The basic problem we wish to consider is- Assuming that X(I) is log-no rmally distributed. determine the probabilit y density functi on I+ en(ln(2z.=., )- rln:'+ (r - l)lnµ ,,)}dx. { of Z(t). Thus we seek to relate the statistics of time.j2(J...:-7Jin'a,, (S) averaged concentrations to those of the raw data In summary, F, is the distribution function of the which are used to construct the averages. We carry out the ana lysis for an averaging period twice the random variable Z(I). which is the average of two length of the fundamental averaging period on which log-normally d istributed random variables. We seek the raw data a re based. In practice, averaging periods to determine how close the distributi on of Z(t) apgreater than twice the basic average are used. For proximates a log-normal distribution. Therefore, we instance, daily averages are computed from 1-h aver- compare Fz from (8) with FI , a log-normal distribuages. The sa lient features of the averaging process tion fun ction with the same mean and variance as will. however. be elucidated with a period twice the F, . Comparisons of F2 and FI were made in 15 cases: a,.,= 1.18. 1.5. 2. 4, 10 and,= 0. 0.5. 0.9. (fhe length of the basic period. We assume that the first order density function of range of a,, for typical air pollution data is between I and 4.) Figures 5 and 6 show comparison of F 2 X(t). px(x ;t). is log-normal. and Fi . Qualitatively we note that Fz compares px(x ;t) = I exp{- {ln x - Inµ,,)'} closely with F1, particu larly at large va lues of :. In v'27; , In a,, 2 ln 'a,, (4) the Appendix we discuss in more detail the compariand that the second order density fun ction px(x. t - r: son between F z and Fr We can conclude from the results shown in Figs. y, t) is the following joint log-normal densit y. 5 and 6 that the distribution of the average of two 1 px(x.t - ,;y.t) = [2n xy ln 'a,, . correlated log-normal varia bles approximates a log- ",,1'' JT="?r {lnx - Inµ ,,)' - 2,{lnx - lnµ ,,)( lny- Inµ ,,)} exp + (In r - In µ.,I' , { 2(1-r 2 )ln 'a,, • fl q, ' whereµ,, and a,, are respectively the geometric mean and the standard geometric deviation of X(I). and , is a correlati on parameter. The firs t orde r density function of Z(I) can be determined from the relation 0.01 "' (I • ' .q ,,o.o (5) OJ 001 ( 2: Pz{z;t) =2J P,(x.t-,;2:-x.t)dx. (6) 0001 T he distribution function F z<:;t) is related to Pz{:: 1) by O.Oi 0 Fz<:;t ) = I: OIO tO ~ 90 99 99.9099.99 FREOUE ~CY D~STRI BUT10N FUN CTION . percent P,<~ : t)d,1 . (7) Fig. 6. Comparison o f F, t v F; for a typica lly high value of o,x 149. 0.9 To•,.. - - - - - - - ~minutes I hour 0.8 8hours ----ldoy 0.7 0.6 0.5 0.4 0.2 0.1 I month Fig. 7. Value of the correlation parameter r necessary for the standard geometric deviation of F; to equal o,b plotted as a function of the standard geometric deviation of the concentration distribution for averaging time T 0 . normal distribution. We now return to the observation (2). If we let T, = r and T• = 2r , then we askGiven a,, = a••, can we, by appropriate choice of the correlation parameter r, have the standard geometric deviation of FJ equal a,•. Figure 7 shows that for the values of a,, and r of interest such an r exists. Although the correlation parameter r is not known for actual data, the results of this section indicate that the observation (2) can be explained simply as a consequence of the near log-normality of the shortaveraging time concentrations. Others have also investigated the effects of averaging time on concentration frequency distributions. There does not appear to be any other work which predicts both the form and parameters of time averaged concentration distributions, however, other relationships for measures of the variance in the distribution versus averaging time have been developed. Saltzman (1970) presented an empirical relationship for the standard geometric deviation. while Shoji and Tsukatani (1973) and Larsen and Peterson (1974) developed relationships based on assumed spectral properties of the concentration time series. ANALYSIS OF FREQUENCY DISTRIBUTIONS Up to this point we have shown that several distributions are capable of representing a number of air pollutant concentration frequency distributions. The log-normal distribution, while not, in fact. a perfect representation ·of the data is a convenient one because the parameters of the distribution can be easily determined from a log-probability plot of the data. Assuming log-normality of the basic data. we then showed that time averages formed from the original data would essentially preserve the log-normality of the data. We now come to the crucial question- wh y do the concentrations tend to be approximately lognormally distributed in the first place. This section is devoted to an allempt to propose possible explanations for this observed phenomenon. Because of the approximate universality of pollutant frequency distributions we would expect the principal factors affecting urban concentration frequency distributions also to exhibit universality. That is, the fundamental characteristics of such a factor relative to its influence on urban pollutant concentration must be approximately the same for all pollutants in all cities. Conversely, a factor which is fundamentally different either from city to city or in its effect on different pollutants does not appear to have a major influence on frequency distributions. Certainly there exists no characteristic similarity over all cities or all pollutants among either the spatial distribution or strength of sources. Similarly. wind direction distributions can be argued not to be a major factor in air pollutant frequency distributions. The concentration of a pollutant near the center of a uniform area source will not be sensitive to wind direction. whereas the concentration at a position near a single point source will be very sensitive to wind direction. Thus the impact of wind direction on concentration can range from negligible to significant and certainly varies from location to location. By such reasoning one is led to the conclusion that the two factors most likely to influence air pollutant frequency distributions are wind speed and mixing height. Increases in both factors will lead to a de- crease in concentration. A number of studies have been carried out investigating the correlation between wind speed, mixing depth, and air pollutant concentrations. Schmidt and Velds (1969) calculated a correlation coefficient between yearly average SO, concentrations in Rollerdam and wind speed of -0.97. Marsh and Withers (1969) found significant correlation between 6 h average SO, concentrations in Reading and wind speed but not between vertical turbulence and concentration. Similarly, in analyzing hourly ozone data in New York City. Bruntz er al. (1 974) determined close correlation with wind speed. and were not able to improve the correlatiori by including mixing height . While certainly not proving that wind speed is the sole factor governing. pollutant frequency distributions, these studies do indicate at least that wind speed/cdncentration correlations are as one might expect. Because of the difficulty in measuring. mixing depths relative to wind speeds. fewer studies exist wherein mixing depth/concentration correlations were computed. Undoubtedly. mixing depth does play a role ii' d! 1erm1ning pollutant frequency distributions. 150. 20 and ~ is the molecular d iffusivity of the species in air. Molecular difTusion is normally neglected when (9) is applied to atmospheric species. the result being the so-called advection equation. o 3-4 minules A, 11.el'I, SC, 1966 ·1968 (Coope, and Rutche. 1968 ) -,:: 10 o 3 hou•s}locQ . E " 1 week ci F,onc e. 1968 · 196~ (Beno,.,, t9 74 1 <'C :::: a. ·-· + V · UC = 0 . (10) i'r "' 0 z §: bo1 010 1 10 FRE QUENCY, percen1 .99 Fig. 8. Wind speed distributions in two cities. Data from Lacq shown for two averaging times. From the above qualitative arguments we have arrived at the conclusion that wind speed and mixing depth should be the major factors influencing air pollutant frequency distributions. Following this supposition, we need to be more quantitative. that is to determine in what manner wind speed and mixing depth influence frequency distributions. The remainder of this section is devoted to an attempt to explain observed frequency distributions on a fundamental basis. In doing so. we restrict our attention to non-chemically reacting pollutants. Because the wind speed U is a random variable, the instantaneous concentration C is a random variable. As noted above, let us consider a one-dimensional flow (in the x-dircction) into which a pollutant is steadily emitted at the x = 0 plane at a rate S g cm · ' s • 1 . The wind velocity in the x-direction is taken to be a function of time only. i.e. U, = U(I). This case, although highly simplified, exhibits the basic features of the situation in which there are three velocity components. The instantaneous concentration is described by the one-dimensional form of (IOI. ?.C(t~ + U(t) il_CU . x) = 0 ar subject to C(O, x) = 0 (12) C(t.0)= ~ r>O . (13) U(rl Influence of wind speed on insrantaneous concentrations In this subsection we wish to consider the influence of wind speed on the distribution of instantaneous pollutant concentrations. To begin, we need some notion as to observed frequency distributions of wind speed. Figure 8 shows frequency distributions of wind speed from Lacq . France and Aiken, SC. In both cases we note that the wind speed is approximately log-normally distributed. Why wind speeds seem to be log-normally distributed and even. in fact, if the log-normal is the best distribution to describe wind speeds are questions not central to our purpose here. We seek oPly to ascertain the consequences of this approximate log-normality in determining pollutant frequency distributions. We begin in some sense with the most basic problem, tha t of determining the effect of wind speed variations on the instantaneous concentration of a pollutant released into that wind field . Clearly. attacking this problem on the basis of a three-dimensional urban flow is impossible. Therefore. we need to isola te the key elements of the problem in a much simpler hypothetical situation. which, nevertheless. retains the basic physics. Such a situation is embodied in a one· dimensional flow into which a pollutant is steadily emitted at a plane. The fundamental equation describing the instantaneous concentration C of an inert atmospheric species is the continuity equation . ~~- + V· lJC = !Y\1 1 C (91 clr where U is the instantaneous wi nd velocity vector. (11) i'x The solution of(l 1- 13) is C(t. x) = ~, V(t) ,. (141 [. u,,·"',· Equation (14) relates the concentration at any position x and time r to the source strength S and the wind speed. We recall that U(t) is a random variable. and therefore that C(r,x) is a random variable. Given the pdf of U(r), Pu(u:t), we wish to determine the pdf of C. pc{c:r.x). It is advantageous to assume that there are only a finite number / of possible wind speeds. with the maximum wind speed denoted by u1 . In addition. we assume thal a given wind speed persists for a time ti.1 . If ti.x is the distance a fluid element moves in ti.r corresponding to the slowest wind speed u 1. i.e. u, = tJ.x/tJ.t. and if the wind speeds obey the relations. u; = iu 1 . i = 1.2 . ... . /, then the proba bility of observing concentration c; at time r, = 1161 and posi\ion -'m = mlix. pc{c;: r,.,.I· is given hy x"'2x 1 and1 1 $ l 11 < Y. I p,:{c,:r , .xml = L P,~uj ; l ,)pc{c;:t , . ,.Xm - j l j c: \ 1.2. 3. .. I • ( I 5al - x. < x"' $ 0 or rr- ~-~ d p,{c;: r,.xml ,'= ,·,= 1. 2..... / . (l5bl The firs! term in (I 5al represents lh< probability of a fluid clement ~i th concentration c; being 151. 0 ~ V) 10 . 0 z 3: d.0 1 OJO-- 1- ·- 1()---l.. ..... ~~ 9 0 - gy- 99909999 tREOUENCY, percen t Fig. 9. Assumed wind speed distrihution for one-dimensional advection problem. advected to location Xm while the second term represents the pro bability of 'fresh· emissions reaching xm . As I is increased, this formulati on yields ·an appro ximation to the continuous density Pc{c :1.x). The solution was evaluated for the case I = 30, !J.r = I and u, = 30. The wind speed was assumed Lo be log-normall y distributed in accordance with observation. Its frequency distribution is shown in Fig. 9. The resulting frequency distribution for the concentration at x = 30. 1 = 30, C(30.30). is shown in Fig. 10. The result is a close approximation to a log-normal distribution. In summary. this situation of a steadily emitting source in a one-dimensional flow in which the wind speed is log-normally distributed leads to an instantaneous concentration that is approximately lognorma lly distributed. This example. while clearly highly idealized. does. however. demonstrate rigorously tha t log-normality of wind speed does lead to log-normalit y of instantaneous concentration. Of course. in an atmospheric fl ow o ther phenomena will influence concentration distributions. and so this example docs not necessaril y establish the cause of approximate log-normality in air qualit y data. Influence cf wi nd speed and mixing d«.'pth on mt>un co11centra1ions In the previo us subsection we considered the effect of wind speed distribution on the instantaneous concentration of a species. Monitoring data reflect the instantaneous concentration at a point. This instantaneous concentration is a random quantity because of the turbulent nature of the atmosphere. When air pollutant concentrations are analyzed from a theoretical point of view. only the mean concentration (c) can be pred icted. where C = (c) + c'. (The mean concentration ( c ) is a function of location and time but is theoretically the result of an ensemble average.) Virtually all mathematical models of air pollutant be• h is reaso nable to suppose that the distribution of mean <:oncentrar ion.,; over a time period of the order of several months to one year is not significantly different havior are concerned with the prediction of (c) (Lamb a nd Seinfeld. 1973). There exist a large number of urban air pollution models depending on the dynamic nature (steady vs unstead y). type o f source (pnint . line, area). number of spalial dimensions, meleorological assumptio ns. boundary conditions, etc. In these models the inputs are usually specified as known. or deterministic. qua ntities (such as wind speed and mixing depth). However. in allempting to assess the effect of the variability of these inputs on the predicted mean concentration. one can propose to allow the inputs to assume a distrihution of values and determine the resulting dis1rihut io11 of values of the mean concen tration . In panicular. we are interested . of course. in the effect of the variability of wind speed and mixing depth on the distribution of mean concentrations. • Box model. II is sometimes ass umed that the mean concentration in a local region of an urban area can be represented by the simple box model relationship. k (16) (c)= ~1,· where u is the mean wind speed. h is the mixing depth, and k is an empirical proportionalit y constant. If we allow u and/or h to be random variables. then the distribution of ( c ) can be computed directl y from (16). Kn ox and Lange (1974) employed this approach by assuming a frequency distribution for u from hourly average readings and computing the frequency distribution of (c). In the cases presented. the approxima te log-normalit y of the observation was reproduced. We have investigated the effect of assuming different form s for the wind speed distribution on the concentration dis1ribution using (16). In addition to the log-normal. the Weibull. gamma. uniform. a nd bounded distributions were investigated (fable 4). In all cases we assumed tha t the mean wind speed was the same and for the two-parameter distributions tha t the wind speed variance was also the same. The wind speed densities along with the resulting concentration densities are given in Table 4. In most cases the high concentration end of the distribution could be approximated by a straight line (although not the one z 0 ;:: <l er S/10 ..,Z u z 8 than the distribution of instantaneous concentrations over YOOa,-,=--!--' t -·...,,_~__._.,_.---'-;s,.----,s,;--c'-""'~ 1 99 the same time period. Therefore. this approach should not be wholly inapplicable 10 the basic issue of analyzing ambient mon itoring data. Fig. 10. Calculat'ed concentration distribution correspond• O.Ot O.IO { .J fR[OUEN~~Y. perce: 99 ing to wind sp«d distribu1ion in Fig. 9. 90 99 • 99. 15cJ. Table 4. Relationship between probability density functions for wind speed and concentration when c = k/u Pdrfor wind speed Pdf for concentration p..(u) !',fr) Log-normal ~ 1 (v' 2nuln/i) " exp { - (In u - In •I'} n 'P 21 Weibull 1 au' exp { - au' ' P+I } Gamma 1 P" I/'(• + I I u' ex p (- pu) Uniform . - I Bounded found by assuming the original log-normal distribution for wind speeds). It is reasonable to conclude that with the simple model of 06) a variety of wind speed distributions could actually exist and still the log-normal would be an adequate approximation to the concentration distribution. The explanation. of course. lies in the inverse relationship between concentration and wind speed. as pointed out previously by Benarie (1969. 1971. 1974) and Pollack (1975). Similarly. the inverse relationship between concentration and mixing depth indicates that the same correspondence between mixing depth and concentration is to be expected . Gaussian plume model. For conditions of (I) a continuous point source located at (0.0,z 0 ) and (2) wind speed constant and direction aligned with the x-axis. the ground-level mean concentration (c(x ..v.0) ) can be estimated by the familiar Gaussian plume equation. 'J}. (c(x,y.0) ) = ~ -- exp{ -:· ,~}exr{=na>.a:u 2a>' 2a, (1 7) Let us assume that the wind speed u is distributed according to a log-normal density with parameters µ, and a,. As expected. we then obtain P,,, ((c)) from ( 17) as log-norma l. I P« ,(( c ) ) = - ,= ·· - - - · ._/ 2n ( c ) In a, 2 ln i a, c'(c ) = -c ( K(z)-8(c ) ) . ilx Cz Oz u(c) - (20) where u(:) = u 1 : ' and K(:) = K ,:1 . The solutions of (20)forground-level crosswind line and area sour= are. (Monin and Yaglom. 1971 ; Lebedeff and Hameed. 1975) and <c< '-°l > = i<.~1~ fl)ri;;,( i~';,;, r 0 ,,, mi where p = • - f3 + 2 > 0. q = (• + I)/(• - f3 + 2). and Os p < I. Equations (21 and 22) can be written in the form <c(x.0)> = Au/ If we now assume u, to be log-normally distributed. we find l'<c,((c) ) to ~e given by (18) with u, = a'.'.' andµ,= Aµ;, . Thus. the eddy diffusion model, while more detailed than the Gaussian model. still predicts a log-normal concentration distribution provided that the wind speed is log-normally, dist ributed. Summar_\· - (In ( c) - In 11, )'} x exp { Eddy difJusion mod<'/. Finally. we consider twodimensional. steady-state diffusion as described by the atmospheric diffusion equation. - , (18) where a,= 0 11 and 2 µ = _ s exp{-\ }exp{- 'f }. ' 'ir.G/1:µu 2a>' 2a; (19) In this sectio n we have attempted to shed some light on the fundament al question of why urban air pollutant concent'ration frequency distributions tend to be log-normdl. Starting from information available from statisti ~J correlations between wind speed and mixing depth (primarily wind speed). we investigated 153. how both instantaneous and mean concentrations might be influenced by wind speed and mixing depth (primarily wind speed) variability. In all cases we found that if wind speeds are nearly log-normally distributed then resulting concentrations will be nearly log-normally distributed. In fact , other distributions. namely the gamma and Weibull. are capable of producing nearly log-normal concentrations. This result does not, of course, establish that wind speeds are the primary influence on concentration distributions in the atmosphere. since other effects are most certainly influential. Nevertheless, the results here are convincing of the role of wind. speed in pollutant frequency distributions. It is interesting to note that in the three simple models considered for the mean concentration. the standard geometric deviation is independent of the source strength and the geometric mean varies linearly with it. Thus. if one were using ( 17. 21. or 22) to predict mean concentrations for source emission changes, only the intercept of the concentration distributions plotted on log-probability paper would change. not the slopes. CONCLUSIONS (Edited by H. M. Englund and W. T. l!crry). Academic Press, New York . Bcnaric M . (1974) The use of the relationship between wind velocity and ambient pollutanl concentralion distrihutions for the estimation of average concentrations from gross meteorological data. Paper 5. Pmn-(•cliny., tlii• S_rmposiwn 011 S1a1i.,1irnl A.'if'('Cls ( f Air Quulir_r Du1t1. Environmental Protection Agency Pub. Nu. ,,r 650/4- 74-0JX. Bruntz S. M .. Cleveland W . S .. Kleiner G. and Warner J. L. (1974) The dependence of amhient ozone on solar radiation. wind . temperature. and mi xi ng height . Symposium on A1r,u1spherl<' Diffusion uncl Air Pollution. Sant.I Barbara. CA. American Meteorological Society. Boston. MA. Cooper R. E. and Rusche B. C. (19681 The SRL meteorological program and off-site dost: calculati ons. E. I. DuPont de Nemours and Co .. Sa\"annah Rivi:r Laboratory. Aiken. SC. Rcpt. DP-1163. Knox J. I!. and Lange R. I 19741 Surface air pollutant concentration frequency distributions : Implications for urhan modeling. J. Air Pollut . Cnntml As.\". 24. 48 SJ Lamb R. G . and Seinfeld J. H. (1973) Mathematical modeling of urban air pollution : General theory . Enriron . Sd. Te,·hno/ 7, 253 -261. Larsen R. I. (1971) A mathematical model for relating air quality measurements to air quality standards. Environ- mental Protection Agency Pub. No. AP-89. Larsen S. E. and Petersen E. L. (19741 Statistical description of air p0llution concentration, averaging. time and frequency. Symposium on Atmospheric Diffwiicm and Air Pollution, Santa Barbara. CA. Amer ican Meteorological Air pollutant concentration frequency distributions Society, Boston. MA. Lebedeff S. A. and Hameed S. (19751 Steady stale solution are the result of complex phenomena. The direct preof the semi•empirical diffusion equation for area sources. diction of these distributions does not appear to be (Preprint). possible. Observed data are generally represented as Lynn D. A. ii 974) Fining curves to urban suspended parlog-normal, although other common statistical distriticulate data. Paper I l Pron't·diny.\ <~{ thl' Sympo.'ii1m1 butions are capable of representing the data as well on S1a1istkal Aspects~( Air Quafil_r Daw . Environmental Protection Agency Pub. No. 650;4- 74-0.1X. as or better than the log-normal. The log-normal is Mage D. T. and On W . R. (19751 An improved statistical convenient because the mean and variance can be model for analyzing air pollution concentration data. easily determined from a log-probability plot of the Paper No. 75-51.4. 6Mt/J Mt'('ling ,f tht· Air Pollution data. The fundamental question of why concentration Cnn1ro/ Association. Bos1on. MA . Marsh K. J. and Withers V. R. (19691 An experimental distributions tend to be approximately log-normal study of the dispersion of the emissions from chimneys cannot be answered unequivocally. The persistence of in Reading- Ill. The investigation of dispersion calculalog-normality for all averaging times can be explained tions. Atmospheric Em:ironmrm 3, 281 - 302. if the raw data are themselves log-normally distriMitchell R. L. (1968) Permanence of the log-normal distributed. The log-normality of the raw data, i.e. the inbution. J . Opt. Soc. Am. 58, 1267- 1~72. Monin A. S. and Yaglom A. M. (1971) Statistical Fluid stantaneous concentrations. can be shown to result Mechanics. MIT Press. Cambridge. MA. if wind speed is log-normally distributed. ConvenPollack R. I. (1975) Studies of pollutant concentratio11 fretional models for mean concentrations. such as the quency distributions. Environmental Protection Agency Gaussian plume and eddy diffusion. can also lead to Pub. No. 650/4-75-004. log-normality for concentration distributions if the . Saltzman B. E. (1970) Significance of sampling time in air monitoring. J. Air Pollut. Control Ass. 20, 660-665. wind speeds are log-normal. It is shown how these Schmidt F. H. and Velds C. A. (19691 On the relation models may be used to estimate the shift in the distri- • between changing meteorological circumstances and the bution resulting from source emission level changes. decrease of sulphur dioxide concentration around Rotterdam. Atmospheric Environment 3, 455---460. Acknowledgement ·- This work was supported by National Science Foundation Grant ENG71-02486. Shoji H. and Tsukatani T. (1973) Statistical model or air pollutant cOncenlration and its applica1i on to the air quality stan¼fards. Armo:,pheric Em·irm,mem 7, 485-501. REFERENCES Benarie M. (1969) Calculation of the amount and the nuisance of pollutants emitted from a point source. Atmospheric Environment 3, 467-473. Benarie M . (1971) About the validity or the log-normal distribution of pollutant concentrations. Paper SU-18D. Proceedings of the 2nd International Clean Air Co1111ress APPENDIX. TH ~ DISTRIBl1TIOI" OF THE AVERAGE OF TWO CO RELATED LOG-NORMAL VARIATES Equa1ionS (7 and 8 1 give the distribution function. f 2 . of the average of two correlated log-normal variates. Figures 5 and 6 show typical comparisons of Fz to t 7 . the log-normal distri bution func1ion with the same mean 154. Table A. I. Comparison of the actual difference F z{:)- F~:) with the first co rrection term E( : ) for the di stribution of the average of two log-n o rmal variates a,, µ,, 1.5 JO £(:) 0.0 0.9 4.0 0.01 0.0 0.9 4.32 7.76 10.4 13.9 25.0 307 6.76 10.0 14.9 32.9 j 0..\97 X JO 0.411 X JO · 0. 1)2 X JO-' 0.42) X )0- O 0.440 0.18 1 X 10 · ' 0.271 X JO 0.)05 X JO ' 0.40) X 10 · ' 0.6(11, and variance as F1 . Qualitati ve comparison such as these indicate that Fz is itself nearly log-normal. Following Mit chell's (1968) analysis of the sum of n independent log-normal vari ates. we investig.ate anal yt ica ll y the near log- no rmality of Fz in this Appendi x. The pdf. p1 . of th e average of two correlated log-normal variat es can be expressed as an orth ogona l polyno mial expansion in terms of p~. the pdf of the log-normal density with the same mean and variance as Pz. i.e. pz(:) = b 0 u0 (:)p;(:) + b,u,(:)p;(:) + (A.I) where the u" are the orthogonal polynomials o r order n defined b} f u.(qlu.(q)p~q) dq = g: : : and the h" a re constants Because pz and Pi ha ve identical mean s and va.riances (A.I I becomes -0_4 I) X JO - • - 0.8)6 X JO · ' 0.231 X w ·· -' 0.829 X 10 -j -0. )29 X 10 · ' - 0.930 X 10 -• -0.)()0 X IO ➔ " -0.778 X JO -' 0.255 X io -• -0.345 X 10 · 7 -0.222 X 10 - " - 0.232 X JO - • -0.57) X 10 · • -0.501 X 10 · • 0.5)9 X JO · ' -0.246 X 10 - l(I - 0.271 X 10 · • -0.727 X JO - 0.771 X JO - ' 0.2 1K x 10-" ' -0.JQ4 X 10 · 3 -0.950 X io · ·• 0.791 X JO (J_t,9(, X JO 0. )07 X IU -0.749 X 10 -0_:\49 X 10 0J .1tJ X 10 -0. IKJ x JO - 0 19 .' x JO -0. 111 X 10 -0.4)0 X 1(1' I -0.2 49 X 10 ' 0.5)9 X 10 -0.494 , 10 -0.4K 4 x 10 - 0.79) X 10 ' -0.557 X 10 0..17K x 10 0J3(1 X JO . The resultin g integra ted expression fo r £(:} is where µ 3 1 = third central momen t Pz µ) , = third central moment of pi :x 1 1 p2 =firs! moment about the origi n o r l'z =~_:~~f Pi = exp lln a•.: medi a n o r Pi '· 1 (A. 8) 1 1 N~:I = ~{erf [ - ~ ( "..:...::-_l_n_i'~: - klna•, )]+ 1}. 2 , ·2 In a~: ' µ~: = geometric mean or Pi a ~1 = standard geometric deviat ion o r pi, or in terms o f the parameters of the initia l St..-cond pzl.:1 = p~:) + h 3 u 3(:)p;l :I + h,u 4 (:)p~:I + (A.4) order densit y run cti o n (e quati ons 3 and 4) In this form we see that r z is the long-normal den sit y. µ,, - µt = 1 [ exp (1 n -", ' .. ) - exp (r 1n ·a,J ' ] •' Pt plus a .. correction .. consisting or an infinite series or (~~7· 8 polynomials multipl ying p~ . Ir we truncate (A .4) after the first correction term and (I) if the truncated expansion is , ( exp (In 'a,, ) )' ' 9 a good a pproximation to Pz and (2) ir the co rrectio n term µ, : - µ,,. jexp(ln 2a,_. I + ½exp tr In .!a,) (A. ) is small. then we see analytically that pz is in fact • a~ 1 = exp :in [l /2ex p(ln 2a,_.) + l,'2 ex p(r In 2a,.,)] : 1 2 . nearly log-normal. Similarly we can integrate the trunca ted expansion and fo rm the same conclu sions abo ut the For the 15 cases ment io ned earlier. £ was evlauated and distribution functi o n F2 . Th e expa nsion is then compared to the actual difference between F2 and F~. Th e result s of fo ur, of these cases are ~iven in Table A.I. These (A.51 result s are typjcal in that they show that (I) fo r low va lue!where of a,,. the difference between F 2 and F2 1s small and L is a good es1imate of this difference and (2) for relati vel y high values or a,J the difft:rence be tween f z and 2 is (A .6) still small but tf, underestimates this difference. r