- Email: [email protected]
Prediction of Pollutant Concentration for Air Pollution Control
PREDICTION OF POLLUTANT CONCENTRATION FOR AIR POLLUTION CONTROL Y.Oshima,* N. Mohri,* H. Nose* and K. Nakagawa** *InsWute of Industrial Sdence, University of Tokyo **Fuyo Data Processing & Systems Development, Ltd.
1. I ntr odu ction
of the performance index and the state equation. The details of a prediction method are described in the following.
In order to solve the air pollution problem thoroughly, it is necessary to find the manufacturing processes which prevent emission of the air pollutants. Until such a final goal is reached, the air pollution should be controlled keeping balance with the production activit y . If the production a ctivity is suppressed extremely to avoid the air pollution, civilized living becomes impossible. On the contrary, if emission of the air pollut ants is not restricted, existence of human being itself becomes dangerous. Therefore, well-balanced control of the air pollution is necessar y. The problem of optimum a ir pollution control is approached from the standpoint of the s ystems engineering and the modern control theories are applied to air pollution control. Here control means that distribution of the pollutant concentration a few hours later is predicted based on the information concerning the pollutant concentration and the me teorc1.ogic al condi.tions sent periodically from several monitoring stations spatially distributed in the area and then the command signals transmitted to emission sources for temporary regulation of plant operation are generated by on-line c omputation based on a performance index including social and economic a ctivities. This paper deals with a prediction s yst em as an important subsystem for air pollution control. As the state equation expressing the field of pollutant diffusion, the Fick's partial differential equation including the stochastic noise is adopted. Discrete approximation for nume rical computation i s done by the weighted residual method which is one of the finite element methods. Since the diffusion coefficients included in the state equation depend upon the atmospher i c turbulence, it is rea sonable to regard them as time variant pa rameters. Therefore, the state equation becomes nonlinear. The unknown diffusion coefficients and the distribution of pollutant concentra·tion of the area can be estimated from the noisy measurements of the pollutant concentration and the wind speed at several monitoring stations by the extended Kalman filter which is a second-order nonlinear filter. If the unknown diffusion coefficients are estimated , the state equation is determined. Thus, optimum control of pollutant emission c an be done on the basis
2. Digitization of di f fusion equation by weighted r esiduaZ method 2.1 Di ffusion equat ion Behaviour of pollutant diffusion is expressed by the following Fick's partial differential equation: N
~~ + (U· V)·c = V·(K· Vc) + L
m=l
~ O(x-xm)
. . . . • . . . . . .• (1)
where c: pollutant concentration, U: wind speed vector, x : position vector, V: Hamilton operator (a/ ax, a / ay, ~ / a z), K: diffusion coefficient (K x , Ky ' Kz), Qm: pollutant emission at emission source (xw) and 0 : Dirac delta function. As for the boundary conditions, the following 4 different conditions are considered. 1) Pollutant concentration on the boundary (cQ) is specified, 2) Flux of pollutant passing through the boundary surface with unit normal vector n [Jt = (-KVc + Uc)nJ is specified, 3) Transportation of pollutant by diffusion through the boundary surface [Jd = (-KVc)nJ is specified and 4) Assuming that transportation of pollutant through the boundary surface is proportional to the difference between pollutant concentrations inside and outside the region and letting e be proportional coefficient and 2 be exterior pollutant concentration, the transporta~ tion [Jdn = e(c-c)J is specified.
2.2 Weighted r esiduaZ me t hod Letting L be any differential oper·ator, the differential equation defined in the domain V is considered. . . • . . . . . • . •. (2)
L [
Multiplication by an a rbitrary function W and integration over the dom a in give the folluwing equation
J W·L [
=
0
• . . • . . . . . . .. (3)
On the contrary, if Eq . (3) holds, the function
645
Y. Oshima et al.
646 cj>
Na ca + Nb cb + Nc Cc + Nd cd
is written as n cj>
=
............
L a. f. 1. 1. i=l
. . . . . . . ..
(4) Eq.7 is expressed by
where ai (i=l to n) : unknown parameters. W is also assumed to be expressed by a linear combination of n functions WI, W2 '" Wn n W=
c
(1, x, y, z)
=
m··········
. . . . . . . . . . .. (5)
L b . W. j=l J J
(1, x, y, z)
L[ L a i fil dV i=l
o
n
..... (6)
Eq.(6) means algebraic equations with n unknowns. Solving Eq.(6) and substituting these solutions into Eq.(4), an approximate solution of Eq.(2) is obtained. The functions Wj are called the weighting functions. If Wj = fi, the method becomes the usual Galerkin approximation.
2.3 Shape f unctions When digitization of a partial differential equation is done by the finite element method, the domain to be considered is divided into finite numbers of elements and the approximating functions fi in Eq.(4) are taken as assembly of the shape functions defined below. In the cases of 3-dimensional problems and division of the domain into the tetrahedrons shown in Fig.l, the shape functions are obtained as follows: The vertexes of the tetrahedron a -d are called the nodes. Let ca - cd be the pollutant concentrations at the nodes. The pollutant concentration at any point within the element is approximated by interpolation from ca - cd. In order to obtain the interpolation functions we assume that the pollutant concentration c can be expressed by a linear combination with respect to x,y and z, that is c =.a + Sx + yy + oz
.. ... .....
-l
1 xb Yb zb
.. (11)
[ 1 Xc Yc Zc
1 xd Yd Zd
Since bj are arbitrary constants, the following equation should exist. L[ L a fil dV = 0 i i=l
~
1 xa Ya za
n
L b. j=l J
(10)
From Eqs.(8), (9) and (10) we obtain
where bj (j=l to n) arbitrary constants. Substituting Eqs.(4) and (5) into Eq.(3) gives n
(9)
(7)
Substitution of the coordinates of the nodes Xa - Zd and the pollutant concentrations at the nodes Ca - Cd gives
The interpolation functions Na-Nd in Eq.(9) are the functions of coordinates x,y and z as shown in Eq.(ll) and are called the shape functions.
2.4 Digitization procedure When the diffusion equation is digitized by the finite element method, the tensor representation is used for simplicity. Suffixing means the components of a physical variable with respect to the spatial coordinates. Repeated suffices in the same term imply summation. The suffix appearing just after the comma means partial differentiation in those directions. Using such a notation, Eq.(l) is expressed as c + Ui c'i = (K i c'i)'i + L Qm o (x-xm) ··(12) where C = dC/dt. Multiplying the both sides by the weithting function Wand integrating over the volume of one element give the following relationship by means of the Green's theorem.
f v HcdV + f v WU.C, .dV + fv H'iKic'i dV 1. 1. f
Since a - 0 are obtained by Eq.(8), substitution of these values gives the following rela tionship c = (N a , Nb' Nc ' Nd)
lIK . c . n.dS + L f WQ o (x-~)dV" (13) 1.'1. 1. m v m
Here the weighting function W is expressed by the shape functions No. (a = a,b,c and d) as W= N
a.
H
(14 )
a
W denotes H at the a node. S~bstituting Eq.(14) into Eq.(13) and considering that Wo. is arbitrary, we obtain
+ f v N0. U1.1. .c, .dV + f N, . K.c, .dV V 0. 1.1. 1.
N cdV
fv
a.
=
. . . . . . .. (8)
s
f s Na. K.c,.n . dS + L f 1. 1. 1. m v
~Na.
o (x-xm)dV
...... ..
(15)
Furthermore, c and U are also expressed by the shape functions NS and ' \ ( S or y = a, b, c and d) as c
U
i
Cs
., . . . . . . . . . . .. (16)
Ny Uyi
... , . . . . . . . . .. (17)
NS
= =
where U . : wind speed in i direction at y node. S~bstituting Eqs.(16) and (17) into
Prediction of pollutant concentration Eq.(15), we obtain AaB c B + Ba S Cs + DaS c B E
a
+ L In N <5 (x-x )dV m
1n a
m
.... (IS)
where the matrices are given as
........
A as
IN
dV a NS
B as
IN
N dV'U a NS y yi
DaS
IK.N
(20) (21)
E = IN K. c'i n. dS a a ~ ~
(22)
Regarding integration of the product of the shape functions Na 1 Nb m Nc n NdP (l - p : integer) with respect to the volume, the following formula exists.
I N 1 N m N n NdP dV b
v a
c
I! m! n! p! (1 + m + n + p + 3)!
ing the node a, summation of the converted diffusion fluxes Ea becomes zero since Ea is defined as integration with respect to the boundary line. As for the node on the boundary, for example, c, summation of Ec, that is, ~ Eg gives the diffusion flux through tte line segment db converted into the node c as
(19)
a , i NS , ~. dV
1
647
6V .. , (23)
L E )J )J c
(27)
The same situation exists in the cases of 3-dimensional problems. Cancellation on the boundary surface between two elements is due to the fact that the normal vector of each element has the opposite direction as shown in Fig.3. Thus, assembling the nodal equations reduces the number of the unknown variables in the assembled equation since summation of the converted diffusion fluxes become zero except the boundaries. Therefore, coupling the boundary conditions with the assembled equation, we can obtain the solution.
3. Extended Kalman filter This formula is applied to integration in Eqs.(19) - (21). Ea is interpreted as the diffusion flux through each surface of a tetrahedron converted into the node a using the weight Na and is called the converted diffusion flux. Differentiation with respect to time in Eq. (IS) is approximated by difference as Cs
{c S
Cs
8c
S
+
c ( 0 ) }/.6t S ( 0) (1- 8 ) Cs
........
(24)
. .......
(25)
where c S ( O) and Cs are the pollutant concentrations at time t and t + 6 t respectively and 8 is a parameter such as 0 ~ 8 ~ 1 and normally 8 ~ 0.7. Substituting Eqs.(24) and (25) into Eq.(lS) gives 1
[ 6 t AaS + 8 (B aS + Das ) } Cs
(6~
AaS + ( 8 - l)(B
aS
+ DaS )]c ( 0 )
+ Ea + m L I v QmNa <5 (x-x m) dV
s
·······(26)
Eq.(26) is the equation with respect to the node a of one element and is called the nodal equation. The equation holds at each vertex (a - d) so that 4 equations exist for one tetrahedron. The unknown variables in Eq. (26) are Cs and Ea (a or B = a-d), whose number is S. Assembling Eq.(26) over the domain and considering the boundar y conditions, the number of equations coincides with the number of unknown variables.
As described above, the state equation (IS) involves the unknown parameter Ki' Therefore, Ki is considered as the state variable together with the pollutant concentration c. Thus, the state equation becomes nonlinear with respect to the extended state variables. Then, estimation of the state variables is done by means of the nonlinear filter suggested by M. Athan. The state equation and the observation equation are written as
z=
fez) + u
..........
(2S)
..........
(29)
where z: state vector, f: nonlinear function of z, Yk: observed output vector at time t=tk, H: observation mantrix, u: Gaussian system noise with zero average and covariance Q and vk: observation noise with zero average and covariance R. The state estimate w between the observation times is given by
+
wet) = f(w(t»
~
n
L i=l
~i
trace [Fi(w(t»S] ........ (30)
with the initial condition (31)
In Eq.(30) vectors
~i
denotes the natural basis
2.5 Assemb led equation Assembling all the nodal equations, the assembled equation governing the domain is obtained. For simplicity the case of a 2dimensional problem consisting of the triangles as shown in Fig.2 is considered. Taking notice of the node a and assembling the nodal equations of the elements involv-
,
~2
6
, ~n 6
rl···
(32)
S is the covariance matrix of the estimate error and Fi is the Hessian nxn symmetric
Y. Oshima et al.
648 matrix whose elements involve
F.
a,S
laS
where ~ is the matrix satisfying the equation
=
A(w(t»S(t) + S(t)AT(w(t»
(42)
In Eq.(40) the difference approximation given by Eq.(24) is adopted and the values at the former step are used for cS(O) and the nonlinear term
S is obtained by solving the following equation Set)
...........
1,2, ...... , n···(33)
+ Q(t)
..... (34) with the initial condition
Thus, the following equations are obtained.
.............
(35)
1
(lit
Aa S + Ba S + Da S)c S
where A is the nxn Jacobian matrix with elements
=,
1
+ - A
A ..
i,j
lJ
1,2,
n
The new covariance matrix Lk+l is obtained from ..... (39)
The digitized diffusion equation (18) was the nodal equation. The assembled equation is of the same form with Eq.(18), but the suffix a denotes serial numbers of all the nodes of the domain to be considered. The state equation regarding the diffusion coefficients
IS.
=
0
is added to Eq.(18). Then, the extended state equations are given by
+
as
L ~ l. trace [A
-1
as
F.S]···· (43) 1
K = K( O)
The new state estimate Zk+l is obtained from
Cs
2
••• (36)
Renewal of the estimate at time tk+l when new observation is made is done as follows: The weighting matrix Gk+l is obtained from
AaS
1 A c ( 0).+ E ( ) aSS a + ~ ~Ct xm
lit
=
1
As for
~
(44)
the difference approximation (45)
is used. Using these equations, the state estimate between successive observations is done. Since the element of the observation matrix Hij is Nj(Xi, Yi, zi), renewal of the state estimate at observation time tk+l is done using this Hij and Eqs.(37) - (39).
4. Programming techniques When the number of the nodes resulted from division of the domain is n, the matrices in Eq.(43) such as Aa S have nxn elements. Therefore, it is necessary to solve the linear equation with n unknowns. If n~500, difficulty arises due to the computing time and the memory capacity. Then, the suitable programming techniques become necessary. As for Eq.(43) the wave front method is used because the matrices such as AaS are sparse. This method was developed for the symmetric matrices by B.M.Irons. In this paper the method is applied to the asymmetric matrices. Let us explain the method in the case of a 2-dimensional problem consisting of the triangles as shown in Fig.4 for simplicity. At first, the following equation is considered.
-(BaS + Da S)c S + Ea + ~Na(xm)
2 AaS
-1
L~i trace [Aa S FiS]
. . . . .. (46) (40) From this equation we obtain the following relationships.
Although the covariance matrix S is governed by Eq.(34), it is obtained from the following equation
S
(47)
649
Pre diction of pollutant concentration Eq.(48) means elimination of Xl. Such a procedure is called the frontal forward elimination. The procedure to obtain Xl by substituting X2 obtained from Eq.(48) into Eq.(47) is called the backward substitution. In the case of solution by the finite element method each nodal equation is assembled using the frontal forward elimination. In such a process F22 + F22 where F;2 is the contribution from the elements not yet assembled should be considered instead of F22 in Eq. (48). Similarly b; should be considered with respect to b 2 . F ll , Fl 2 and F21 other than F22 in Fq.(46) are the contributions from the elements already assembled. In this case adoption of F22 + F; 2 instead of F22 in Eq.(46) is equivalent to addition of F; 2 to the lefthand side matrix in Eq. (48). This means that the frontal forward elimination needs not the complete assembled matrix. In the case as shown in Fig.4 the wave front method is done as follows: (1) The element is governed by the following equation.
CD
tt :t: :\j [:;] [:\]
(49)
Since there is no variable to be eliminated, computation proceeds to the element (2)
l:h" :~dal 7"[i~]Of
~;:
GD t[h: 2]element GD
:1: :;: :::1
(50 )
Since Xl does not appear hereafter, it can be eliminated. In this case, however, it is necessary to store the data required for the backward substitution such as -F~: Fl 2 and F~i b l in Eq.(47) in the working area. Thus, the 1st row and the 1st column in Eq. (51) become unnecessary and the following equation results.
a4 22 4 a 62 a4 52
IFA C S.E.S.- W·
a4 26 4 a 56 a4 56
a4 25 a4 65 a4 55
0x
2
x x
5 5
0b 4
2 b4 6 b4 5
....
s
as 55 55 as as 55 55 as as 105 105
t
1010
la
is
(53)
la
When assembling Eqs.(52) and (53) the area for the new variable XIO is allocated to the empty area in Eq.(52). Thus, the following equation is obtained b5 X a5 a5 a5 a5 la la 105 105 1010 102 5 5 5 b 5 5 X a a a a 2 2 25 ··(54) 25 22 210 5 5 5 5 b 5 X a a a a 5 5 55 55 52 510 5 5 5 5 b 5 X a a a a 5 5 55 55 52 510 (4) Similarly, when the nodal equation of the element is assembled, Xs and x9 are eliminated resulting in
®
a9 105 a9 25 a9 55
9 a9 1010 102 9 a9 22 210 a9 52
X
b9
x
b9 2 b9 5
la
x
2 5
0
0
la ··(55)
0 0
The size of the frontal matrix in Eq.(55) is 3. In such a manner the elimination process is repeated until the last element is reached. After that, the backward s ubstitution is done using the data stored in the working area and in consequence all the solutions are obtained. In the case of the example as shown in Fig. 4 the maximum size of the frontal matrix is 5 and the forward elimination is completed by the area 5x5. When dealing with, for example, 500 nodes, the area for the frontal matrix uses about 2.5KW storage capacity and the working area uses about 25KW storage capacity. On the contrary, since the matrices ne~ess ary for filtering are not sparse, the above-mentioned method can not be applied. Therefore, the number of the nodes where filtering ~s applied is limited to about 100. Hence, the nodes of the domain are divided into the nodes where filtering is applied or not as shown in Fig.5. At first, the nodal equations of the elements to which filtering is not applied [The elements with * mark in Fig.5] are assembled as
(~ A' (52)
0)
:i::]n- [:;]. . .
is
When Eqs. (49) and (50) are assembled, the areas with respect to Xl and x 6 are already reserved in the frontal matrix. It is only necessary to expand the are a for xs· Hence, the following equation is obtained. b3 X a3 a3 a3 1 1 15 12 15 11 3 3 3 b3 x a a a 2 2 25 22 25 21 .. (51) b3 X a3 a3 a3 a3 5 55 5 52 55 51 X b3 a3 a3 a3 a3 5 5 55 52 55 51
I
where (====:J means the empty area. (3) The nodal equation of the element
lit
as
+
as + D~S ) cS
B'
'" ~ A' Cs( 0) + E' + li t as a
~
Na (xa )
...
(56)
Y. Oshima et al .
650
The pollutant concentrations Cs involve those at the nodes with 0 mark and those at the nodes with • mark which are on the boundaries between the elements with * mark and the elements with hatching. Eq.(56) is expressed as
[::1- ~:J [::1
. . .. (57)
+
where cl is the concentrations at the nodes 'vi th 0 mark, c 2 is those a t the nodes wi th • mark and E2 consists of E~ with respect to the nodes on the boundaries between the filtered and unfiltered elements. Rll - R22 are the matrices depending on A~S' B~ S and D~ S and El' Tl a nd T2 are given. Elimination of Cl from Eq .(57) gives -1
(R 22 - R2 1R11 R12)C2 -1
E2 + T2 - R21Rll (El + T1) . . ....
(58)
Eq.(58) is the equation governing the nod es on the boundaries between the filtered and unfiltered elements. Then, th e assembled equation with respect to all the filtered e lements with hatching is forme d and the equation number (18) is given to thi s assembled equation. Summation of Ea in Eq.(18) and the corresponding term of E2 in Eq.(57) results in zero by cancellation as mentioned in Sec. 2.5. Hence, Ea in Eq.(18) can be eliminated using Eq.(57). The concentration s at the nodes with 0 mark are obtained by substitution into Eq.(57).
5. Computational results A prediction method described above was applied to the case study dealing with S02 pollution problem in the medium siz e city Himeji. The area covering the whole city stretches 23 Km from north to south and 28 Km from east to west. The height to be considered was assumed to be about 6 Km. The area was divided into the meshes as shown in Fig.6. In the vertical direction division into 3 layers whose lowest one is 340 m high, middle one is 690 m high and highest one is 5800 m high was done. Thus, the domain was divided into the hexahedrons. One hexahedron was divid ed into 5 tetrahedrons. In this paper division of the hexahedron was made in 2 ways as shown in Fig.7 and the concentrations at the nodes were computed from averaging both values obtained by 2 kinds of division. Digitization was done by division into 1560 elements and 476 nodes. The time step of the difference approximation was 20 min. The pollutant emission from the main 10 factories indicated by • mark as ShO~l in Fig.6 was allotted to the appropriate nodes according to their positions and their chimney heights. The accurate information regarding the emission strength from hour to hour was hard to obtain, so that the average valu e obtain e d f rom Clle annual total emission quantity was used. There were 9
monitoring stations as shown in Fig.6. The observed data obtained on Jan. 8, 1973 were used. The computational results as shown in Fig. 8 give the distributions of the predicted pollutant concentrations at 5 and 6 oclock A.~1. after renewal of the state estimate made every hour from start at 12 oclock midnight and the modified distributions after renewal by the latest observed data. For this computation the HiTAC 8800 computer was used and the computing time was 3 min and 40 sec. Fig.9 shows the observed concentration data at the monitoring station and the time-dependent variation of the predicted concentrations over a period of 20 hours at the nodes 89, 90 and 98 surrounding the monitoring s tation as shown in Fig. 6. In this case the computing time was 18 min and 37 sec. It should be noted that the moving average is taken for the computational re sults . Since it is not allowed that th e pollutant moves over one element at one time step of the difference, th ere exists the upper limit of the wind speed to be applied according to the t i me step a nd the size of th e e lement. Genera lly speaking, the lower the ave rage wind speed is, the higher the pollutant concentration is. Therefore, the computationa l results as shown in Figs.(8) and (9) deal with the case where the wind speed is negligibly small.
CD
CD
6. ConcZusions In this paper a method of prediction of the pollutant concentration for air pollution control is presented, in which the finite e lement method is used for digitization of the diffusion e~uation a nd the extended Kalman filter is used for es timat ion of the diffusion coefficients and the pollutant concen tration. Further, the programming techniques and several computational results are presented. ,\ n "l'l'r O'1c l1 to tile opt i mum control of air pollution is difficult a t present because the emission strength data of the main emiss ion sources from hour to hour are very hard to obtain. REFERENCES (1) C.C.Shieh, A generalized urban air pollutic n model and its application in the ST. Louis metropolitan area,
IBM Research (1973). (2) J.Kondo, Atmospheric pollution by emission from automobile engines,
Theoretical and Applied 11echanics, 21(1973) . (3) O.C.Zienkiewicz and C.J.Parekh, Transient fie ld problems; Two-dimensional analysis by isoparametric finite elements,
International JournaZ for Nume r ical Method in Engineer ing 2, 61 (1970). (4) Zienkiewicz, O.C. (1971)
The Finite !;'lement Method in Engineering Science, McGraw-Hill, New York.
Prediction of pollutant concentration
651
(5) Whiteman, J.R. (1973)
The Mathematias of Finite EZements and AppZiaations, Academic Press, New York. (6) R.E.Kalman, A new approach to linear filtering and prediction problems, T~ns. ASME, J. of Basia Eng., 82, 35 (1960). (7) M.Athans, R.P.Wishner and A.Bertoni, Suboptimal state estimation for continuous time non1inear system from discrete noisy measurements, IEEE Trans., AC-13, 504 (1968).
d ~----~r-----~r-----~b
f~------~------~~------~c
a
~------~------~~----~d
Fig.2. Converted diffusion flux on boundary line
b
Fig.l. Tetrahedron element
d'
e e a
bL------~
bL-----~
c
Fig.3. Converted diffusion flux on boundary surface
r -_ _ _ _ _ _~------~------~12
*
*
*
*
*
*
*
*
*
*
*
*
~----~~----~~------~8
Fig.4. Two-dimensional problem
*
*
• Filtered nodes
?: Filtered elements
*
*
0 Unfil tered nodes * Unfiltered elements
Fig.5. Application of filtering
652
Y. Oshima e t a l.
2 (1.2.4.5) (2.3.4.7) (2.5.6.7) (4.5 . 7.8 ) (2.4.5. 7 )
2
0.2.3. 6) 0 . 5. 6. 8 ) (3.6.7.8) 0.3.4.8 ) 0 . 3.6.8)
Fig.7. Divisio n of hexahe dron
Fig.6. Divisio n of Himeji city area
(a) Predic ted distrib ution at 5 oclock A.M.
(c) Predic ted distrib ution at 6 oclock A.M.
(b) Modifi ed distrib ution at 5 oclock A.M.
(d) Modifi ed distrib ution at 6 oclock A.M.
Fig.8. Predic ted and modifie d distrib utions of polluta nt concen tration
Prediction of pollutant concentration
Initial Conditions: S02 concentration 12 ppb Covariance of concentration 1.0 Diffusion coefficients Kx 200 m2/sec Ky = 200 m2/sec Kz = 280 m2 /sec Covariance of diffusion coefficients 10.0
ppb • Observed data
30
653
gComputational results
20
10
0 + - - - - - _ T - -_ _--~--~_4~_T--~--~--~_4--_T1--~--+_--r_~--_+------+_--~~
o
2
4
6
8
10
12
14
16
Fig.9. Observed concentration data and predic~ed concentrations at nodes surrounding monitoring station
18
20
h
1. I ntr odu ction
of the performance index and the state equation. The details of a prediction method are described in the following.
In order to solve the air pollution problem thoroughly, it is necessary to find the manufacturing processes which prevent emission of the air pollutants. Until such a final goal is reached, the air pollution should be controlled keeping balance with the production activit y . If the production a ctivity is suppressed extremely to avoid the air pollution, civilized living becomes impossible. On the contrary, if emission of the air pollut ants is not restricted, existence of human being itself becomes dangerous. Therefore, well-balanced control of the air pollution is necessar y. The problem of optimum a ir pollution control is approached from the standpoint of the s ystems engineering and the modern control theories are applied to air pollution control. Here control means that distribution of the pollutant concentration a few hours later is predicted based on the information concerning the pollutant concentration and the me teorc1.ogic al condi.tions sent periodically from several monitoring stations spatially distributed in the area and then the command signals transmitted to emission sources for temporary regulation of plant operation are generated by on-line c omputation based on a performance index including social and economic a ctivities. This paper deals with a prediction s yst em as an important subsystem for air pollution control. As the state equation expressing the field of pollutant diffusion, the Fick's partial differential equation including the stochastic noise is adopted. Discrete approximation for nume rical computation i s done by the weighted residual method which is one of the finite element methods. Since the diffusion coefficients included in the state equation depend upon the atmospher i c turbulence, it is rea sonable to regard them as time variant pa rameters. Therefore, the state equation becomes nonlinear. The unknown diffusion coefficients and the distribution of pollutant concentra·tion of the area can be estimated from the noisy measurements of the pollutant concentration and the wind speed at several monitoring stations by the extended Kalman filter which is a second-order nonlinear filter. If the unknown diffusion coefficients are estimated , the state equation is determined. Thus, optimum control of pollutant emission c an be done on the basis
2. Digitization of di f fusion equation by weighted r esiduaZ method 2.1 Di ffusion equat ion Behaviour of pollutant diffusion is expressed by the following Fick's partial differential equation: N
~~ + (U· V)·c = V·(K· Vc) + L
m=l
~ O(x-xm)
. . . . • . . . . . .• (1)
where c: pollutant concentration, U: wind speed vector, x : position vector, V: Hamilton operator (a/ ax, a / ay, ~ / a z), K: diffusion coefficient (K x , Ky ' Kz), Qm: pollutant emission at emission source (xw) and 0 : Dirac delta function. As for the boundary conditions, the following 4 different conditions are considered. 1) Pollutant concentration on the boundary (cQ) is specified, 2) Flux of pollutant passing through the boundary surface with unit normal vector n [Jt = (-KVc + Uc)nJ is specified, 3) Transportation of pollutant by diffusion through the boundary surface [Jd = (-KVc)nJ is specified and 4) Assuming that transportation of pollutant through the boundary surface is proportional to the difference between pollutant concentrations inside and outside the region and letting e be proportional coefficient and 2 be exterior pollutant concentration, the transporta~ tion [Jdn = e(c-c)J is specified.
2.2 Weighted r esiduaZ me t hod Letting L be any differential oper·ator, the differential equation defined in the domain V is considered. . . • . . . . . • . •. (2)
L [
Multiplication by an a rbitrary function W and integration over the dom a in give the folluwing equation
J W·L [
=
0
• . . • . . . . . . .. (3)
On the contrary, if Eq . (3) holds, the function
645
Y. Oshima et al.
646 cj>
Na ca + Nb cb + Nc Cc + Nd cd
is written as n cj>
=
............
L a. f. 1. 1. i=l
. . . . . . . ..
(4) Eq.7 is expressed by
where ai (i=l to n) : unknown parameters. W is also assumed to be expressed by a linear combination of n functions WI, W2 '" Wn n W=
c
(1, x, y, z)
=
m··········
. . . . . . . . . . .. (5)
L b . W. j=l J J
(1, x, y, z)
L[ L a i fil dV i=l
o
n
..... (6)
Eq.(6) means algebraic equations with n unknowns. Solving Eq.(6) and substituting these solutions into Eq.(4), an approximate solution of Eq.(2) is obtained. The functions Wj are called the weighting functions. If Wj = fi, the method becomes the usual Galerkin approximation.
2.3 Shape f unctions When digitization of a partial differential equation is done by the finite element method, the domain to be considered is divided into finite numbers of elements and the approximating functions fi in Eq.(4) are taken as assembly of the shape functions defined below. In the cases of 3-dimensional problems and division of the domain into the tetrahedrons shown in Fig.l, the shape functions are obtained as follows: The vertexes of the tetrahedron a -d are called the nodes. Let ca - cd be the pollutant concentrations at the nodes. The pollutant concentration at any point within the element is approximated by interpolation from ca - cd. In order to obtain the interpolation functions we assume that the pollutant concentration c can be expressed by a linear combination with respect to x,y and z, that is c =.a + Sx + yy + oz
.. ... .....
-l
1 xb Yb zb
.. (11)
[ 1 Xc Yc Zc
1 xd Yd Zd
Since bj are arbitrary constants, the following equation should exist. L[ L a fil dV = 0 i i=l
~
1 xa Ya za
n
L b. j=l J
(10)
From Eqs.(8), (9) and (10) we obtain
where bj (j=l to n) arbitrary constants. Substituting Eqs.(4) and (5) into Eq.(3) gives n
(9)
(7)
Substitution of the coordinates of the nodes Xa - Zd and the pollutant concentrations at the nodes Ca - Cd gives
The interpolation functions Na-Nd in Eq.(9) are the functions of coordinates x,y and z as shown in Eq.(ll) and are called the shape functions.
2.4 Digitization procedure When the diffusion equation is digitized by the finite element method, the tensor representation is used for simplicity. Suffixing means the components of a physical variable with respect to the spatial coordinates. Repeated suffices in the same term imply summation. The suffix appearing just after the comma means partial differentiation in those directions. Using such a notation, Eq.(l) is expressed as c + Ui c'i = (K i c'i)'i + L Qm o (x-xm) ··(12) where C = dC/dt. Multiplying the both sides by the weithting function Wand integrating over the volume of one element give the following relationship by means of the Green's theorem.
f v HcdV + f v WU.C, .dV + fv H'iKic'i dV 1. 1. f
Since a - 0 are obtained by Eq.(8), substitution of these values gives the following rela tionship c = (N a , Nb' Nc ' Nd)
lIK . c . n.dS + L f WQ o (x-~)dV" (13) 1.'1. 1. m v m
Here the weighting function W is expressed by the shape functions No. (a = a,b,c and d) as W= N
a.
H
(14 )
a
W denotes H at the a node. S~bstituting Eq.(14) into Eq.(13) and considering that Wo. is arbitrary, we obtain
+ f v N0. U1.1. .c, .dV + f N, . K.c, .dV V 0. 1.1. 1.
N cdV
fv
a.
=
. . . . . . .. (8)
s
f s Na. K.c,.n . dS + L f 1. 1. 1. m v
~Na.
o (x-xm)dV
...... ..
(15)
Furthermore, c and U are also expressed by the shape functions NS and ' \ ( S or y = a, b, c and d) as c
U
i
Cs
., . . . . . . . . . . .. (16)
Ny Uyi
... , . . . . . . . . .. (17)
NS
= =
where U . : wind speed in i direction at y node. S~bstituting Eqs.(16) and (17) into
Prediction of pollutant concentration Eq.(15), we obtain AaB c B + Ba S Cs + DaS c B E
a
+ L In N <5 (x-x )dV m
1n a
m
.... (IS)
where the matrices are given as
........
A as
IN
dV a NS
B as
IN
N dV'U a NS y yi
DaS
IK.N
(20) (21)
E = IN K. c'i n. dS a a ~ ~
(22)
Regarding integration of the product of the shape functions Na 1 Nb m Nc n NdP (l - p : integer) with respect to the volume, the following formula exists.
I N 1 N m N n NdP dV b
v a
c
I! m! n! p! (1 + m + n + p + 3)!
ing the node a, summation of the converted diffusion fluxes Ea becomes zero since Ea is defined as integration with respect to the boundary line. As for the node on the boundary, for example, c, summation of Ec, that is, ~ Eg gives the diffusion flux through tte line segment db converted into the node c as
(19)
a , i NS , ~. dV
1
647
6V .. , (23)
L E )J )J c
(27)
The same situation exists in the cases of 3-dimensional problems. Cancellation on the boundary surface between two elements is due to the fact that the normal vector of each element has the opposite direction as shown in Fig.3. Thus, assembling the nodal equations reduces the number of the unknown variables in the assembled equation since summation of the converted diffusion fluxes become zero except the boundaries. Therefore, coupling the boundary conditions with the assembled equation, we can obtain the solution.
3. Extended Kalman filter This formula is applied to integration in Eqs.(19) - (21). Ea is interpreted as the diffusion flux through each surface of a tetrahedron converted into the node a using the weight Na and is called the converted diffusion flux. Differentiation with respect to time in Eq. (IS) is approximated by difference as Cs
{c S
Cs
8c
S
+
c ( 0 ) }/.6t S ( 0) (1- 8 ) Cs
........
(24)
. .......
(25)
where c S ( O) and Cs are the pollutant concentrations at time t and t + 6 t respectively and 8 is a parameter such as 0 ~ 8 ~ 1 and normally 8 ~ 0.7. Substituting Eqs.(24) and (25) into Eq.(lS) gives 1
[ 6 t AaS + 8 (B aS + Das ) } Cs
(6~
AaS + ( 8 - l)(B
aS
+ DaS )]c ( 0 )
+ Ea + m L I v QmNa <5 (x-x m) dV
s
·······(26)
Eq.(26) is the equation with respect to the node a of one element and is called the nodal equation. The equation holds at each vertex (a - d) so that 4 equations exist for one tetrahedron. The unknown variables in Eq. (26) are Cs and Ea (a or B = a-d), whose number is S. Assembling Eq.(26) over the domain and considering the boundar y conditions, the number of equations coincides with the number of unknown variables.
As described above, the state equation (IS) involves the unknown parameter Ki' Therefore, Ki is considered as the state variable together with the pollutant concentration c. Thus, the state equation becomes nonlinear with respect to the extended state variables. Then, estimation of the state variables is done by means of the nonlinear filter suggested by M. Athan. The state equation and the observation equation are written as
z=
fez) + u
..........
(2S)
..........
(29)
where z: state vector, f: nonlinear function of z, Yk: observed output vector at time t=tk, H: observation mantrix, u: Gaussian system noise with zero average and covariance Q and vk: observation noise with zero average and covariance R. The state estimate w between the observation times is given by
+
wet) = f(w(t»
~
n
L i=l
~i
trace [Fi(w(t»S] ........ (30)
with the initial condition (31)
In Eq.(30) vectors
~i
denotes the natural basis
2.5 Assemb led equation Assembling all the nodal equations, the assembled equation governing the domain is obtained. For simplicity the case of a 2dimensional problem consisting of the triangles as shown in Fig.2 is considered. Taking notice of the node a and assembling the nodal equations of the elements involv-
,
~2
6
, ~n 6
rl···
(32)
S is the covariance matrix of the estimate error and Fi is the Hessian nxn symmetric
Y. Oshima et al.
648 matrix whose elements involve
F.
a,S
laS
where ~ is the matrix satisfying the equation
=
A(w(t»S(t) + S(t)AT(w(t»
(42)
In Eq.(40) the difference approximation given by Eq.(24) is adopted and the values at the former step are used for cS(O) and the nonlinear term
S is obtained by solving the following equation Set)
...........
1,2, ...... , n···(33)
+ Q(t)
..... (34) with the initial condition
Thus, the following equations are obtained.
.............
(35)
1
(lit
Aa S + Ba S + Da S)c S
where A is the nxn Jacobian matrix with elements
=,
1
+ - A
A ..
i,j
lJ
1,2,
n
The new covariance matrix Lk+l is obtained from ..... (39)
The digitized diffusion equation (18) was the nodal equation. The assembled equation is of the same form with Eq.(18), but the suffix a denotes serial numbers of all the nodes of the domain to be considered. The state equation regarding the diffusion coefficients
IS.
=
0
is added to Eq.(18). Then, the extended state equations are given by
+
as
L ~ l. trace [A
-1
as
F.S]···· (43) 1
K = K( O)
The new state estimate Zk+l is obtained from
Cs
2
••• (36)
Renewal of the estimate at time tk+l when new observation is made is done as follows: The weighting matrix Gk+l is obtained from
AaS
1 A c ( 0).+ E ( ) aSS a + ~ ~Ct xm
lit
=
1
As for
~
(44)
the difference approximation (45)
is used. Using these equations, the state estimate between successive observations is done. Since the element of the observation matrix Hij is Nj(Xi, Yi, zi), renewal of the state estimate at observation time tk+l is done using this Hij and Eqs.(37) - (39).
4. Programming techniques When the number of the nodes resulted from division of the domain is n, the matrices in Eq.(43) such as Aa S have nxn elements. Therefore, it is necessary to solve the linear equation with n unknowns. If n~500, difficulty arises due to the computing time and the memory capacity. Then, the suitable programming techniques become necessary. As for Eq.(43) the wave front method is used because the matrices such as AaS are sparse. This method was developed for the symmetric matrices by B.M.Irons. In this paper the method is applied to the asymmetric matrices. Let us explain the method in the case of a 2-dimensional problem consisting of the triangles as shown in Fig.4 for simplicity. At first, the following equation is considered.
-(BaS + Da S)c S + Ea + ~Na(xm)
2 AaS
-1
L~i trace [Aa S FiS]
. . . . .. (46) (40) From this equation we obtain the following relationships.
Although the covariance matrix S is governed by Eq.(34), it is obtained from the following equation
S
(47)
649
Pre diction of pollutant concentration Eq.(48) means elimination of Xl. Such a procedure is called the frontal forward elimination. The procedure to obtain Xl by substituting X2 obtained from Eq.(48) into Eq.(47) is called the backward substitution. In the case of solution by the finite element method each nodal equation is assembled using the frontal forward elimination. In such a process F22 + F22 where F;2 is the contribution from the elements not yet assembled should be considered instead of F22 in Eq. (48). Similarly b; should be considered with respect to b 2 . F ll , Fl 2 and F21 other than F22 in Fq.(46) are the contributions from the elements already assembled. In this case adoption of F22 + F; 2 instead of F22 in Eq.(46) is equivalent to addition of F; 2 to the lefthand side matrix in Eq. (48). This means that the frontal forward elimination needs not the complete assembled matrix. In the case as shown in Fig.4 the wave front method is done as follows: (1) The element is governed by the following equation.
CD
tt :t: :\j [:;] [:\]
(49)
Since there is no variable to be eliminated, computation proceeds to the element (2)
l:h" :~dal 7"[i~]Of
~;:
GD t[h: 2]element GD
:1: :;: :::1
(50 )
Since Xl does not appear hereafter, it can be eliminated. In this case, however, it is necessary to store the data required for the backward substitution such as -F~: Fl 2 and F~i b l in Eq.(47) in the working area. Thus, the 1st row and the 1st column in Eq. (51) become unnecessary and the following equation results.
a4 22 4 a 62 a4 52
IFA C S.E.S.- W·
a4 26 4 a 56 a4 56
a4 25 a4 65 a4 55
0x
2
x x
5 5
0b 4
2 b4 6 b4 5
....
s
as 55 55 as as 55 55 as as 105 105
t
1010
la
is
(53)
la
When assembling Eqs.(52) and (53) the area for the new variable XIO is allocated to the empty area in Eq.(52). Thus, the following equation is obtained b5 X a5 a5 a5 a5 la la 105 105 1010 102 5 5 5 b 5 5 X a a a a 2 2 25 ··(54) 25 22 210 5 5 5 5 b 5 X a a a a 5 5 55 55 52 510 5 5 5 5 b 5 X a a a a 5 5 55 55 52 510 (4) Similarly, when the nodal equation of the element is assembled, Xs and x9 are eliminated resulting in
®
a9 105 a9 25 a9 55
9 a9 1010 102 9 a9 22 210 a9 52
X
b9
x
b9 2 b9 5
la
x
2 5
0
0
la ··(55)
0 0
The size of the frontal matrix in Eq.(55) is 3. In such a manner the elimination process is repeated until the last element is reached. After that, the backward s ubstitution is done using the data stored in the working area and in consequence all the solutions are obtained. In the case of the example as shown in Fig. 4 the maximum size of the frontal matrix is 5 and the forward elimination is completed by the area 5x5. When dealing with, for example, 500 nodes, the area for the frontal matrix uses about 2.5KW storage capacity and the working area uses about 25KW storage capacity. On the contrary, since the matrices ne~ess ary for filtering are not sparse, the above-mentioned method can not be applied. Therefore, the number of the nodes where filtering ~s applied is limited to about 100. Hence, the nodes of the domain are divided into the nodes where filtering is applied or not as shown in Fig.5. At first, the nodal equations of the elements to which filtering is not applied [The elements with * mark in Fig.5] are assembled as
(~ A' (52)
0)
:i::]n- [:;]. . .
is
When Eqs. (49) and (50) are assembled, the areas with respect to Xl and x 6 are already reserved in the frontal matrix. It is only necessary to expand the are a for xs· Hence, the following equation is obtained. b3 X a3 a3 a3 1 1 15 12 15 11 3 3 3 b3 x a a a 2 2 25 22 25 21 .. (51) b3 X a3 a3 a3 a3 5 55 5 52 55 51 X b3 a3 a3 a3 a3 5 5 55 52 55 51
I
where (====:J means the empty area. (3) The nodal equation of the element
lit
as
+
as + D~S ) cS
B'
'" ~ A' Cs( 0) + E' + li t as a
~
Na (xa )
...
(56)
Y. Oshima et al .
650
The pollutant concentrations Cs involve those at the nodes with 0 mark and those at the nodes with • mark which are on the boundaries between the elements with * mark and the elements with hatching. Eq.(56) is expressed as
[::1- ~:J [::1
. . .. (57)
+
where cl is the concentrations at the nodes 'vi th 0 mark, c 2 is those a t the nodes wi th • mark and E2 consists of E~ with respect to the nodes on the boundaries between the filtered and unfiltered elements. Rll - R22 are the matrices depending on A~S' B~ S and D~ S and El' Tl a nd T2 are given. Elimination of Cl from Eq .(57) gives -1
(R 22 - R2 1R11 R12)C2 -1
E2 + T2 - R21Rll (El + T1) . . ....
(58)
Eq.(58) is the equation governing the nod es on the boundaries between the filtered and unfiltered elements. Then, th e assembled equation with respect to all the filtered e lements with hatching is forme d and the equation number (18) is given to thi s assembled equation. Summation of Ea in Eq.(18) and the corresponding term of E2 in Eq.(57) results in zero by cancellation as mentioned in Sec. 2.5. Hence, Ea in Eq.(18) can be eliminated using Eq.(57). The concentration s at the nodes with 0 mark are obtained by substitution into Eq.(57).
5. Computational results A prediction method described above was applied to the case study dealing with S02 pollution problem in the medium siz e city Himeji. The area covering the whole city stretches 23 Km from north to south and 28 Km from east to west. The height to be considered was assumed to be about 6 Km. The area was divided into the meshes as shown in Fig.6. In the vertical direction division into 3 layers whose lowest one is 340 m high, middle one is 690 m high and highest one is 5800 m high was done. Thus, the domain was divided into the hexahedrons. One hexahedron was divid ed into 5 tetrahedrons. In this paper division of the hexahedron was made in 2 ways as shown in Fig.7 and the concentrations at the nodes were computed from averaging both values obtained by 2 kinds of division. Digitization was done by division into 1560 elements and 476 nodes. The time step of the difference approximation was 20 min. The pollutant emission from the main 10 factories indicated by • mark as ShO~l in Fig.6 was allotted to the appropriate nodes according to their positions and their chimney heights. The accurate information regarding the emission strength from hour to hour was hard to obtain, so that the average valu e obtain e d f rom Clle annual total emission quantity was used. There were 9
monitoring stations as shown in Fig.6. The observed data obtained on Jan. 8, 1973 were used. The computational results as shown in Fig. 8 give the distributions of the predicted pollutant concentrations at 5 and 6 oclock A.~1. after renewal of the state estimate made every hour from start at 12 oclock midnight and the modified distributions after renewal by the latest observed data. For this computation the HiTAC 8800 computer was used and the computing time was 3 min and 40 sec. Fig.9 shows the observed concentration data at the monitoring station and the time-dependent variation of the predicted concentrations over a period of 20 hours at the nodes 89, 90 and 98 surrounding the monitoring s tation as shown in Fig. 6. In this case the computing time was 18 min and 37 sec. It should be noted that the moving average is taken for the computational re sults . Since it is not allowed that th e pollutant moves over one element at one time step of the difference, th ere exists the upper limit of the wind speed to be applied according to the t i me step a nd the size of th e e lement. Genera lly speaking, the lower the ave rage wind speed is, the higher the pollutant concentration is. Therefore, the computationa l results as shown in Figs.(8) and (9) deal with the case where the wind speed is negligibly small.
CD
CD
6. ConcZusions In this paper a method of prediction of the pollutant concentration for air pollution control is presented, in which the finite e lement method is used for digitization of the diffusion e~uation a nd the extended Kalman filter is used for es timat ion of the diffusion coefficients and the pollutant concen tration. Further, the programming techniques and several computational results are presented. ,\ n "l'l'r O'1c l1 to tile opt i mum control of air pollution is difficult a t present because the emission strength data of the main emiss ion sources from hour to hour are very hard to obtain. REFERENCES (1) C.C.Shieh, A generalized urban air pollutic n model and its application in the ST. Louis metropolitan area,
IBM Research (1973). (2) J.Kondo, Atmospheric pollution by emission from automobile engines,
Theoretical and Applied 11echanics, 21(1973) . (3) O.C.Zienkiewicz and C.J.Parekh, Transient fie ld problems; Two-dimensional analysis by isoparametric finite elements,
International JournaZ for Nume r ical Method in Engineer ing 2, 61 (1970). (4) Zienkiewicz, O.C. (1971)
The Finite !;'lement Method in Engineering Science, McGraw-Hill, New York.
Prediction of pollutant concentration
651
(5) Whiteman, J.R. (1973)
The Mathematias of Finite EZements and AppZiaations, Academic Press, New York. (6) R.E.Kalman, A new approach to linear filtering and prediction problems, T~ns. ASME, J. of Basia Eng., 82, 35 (1960). (7) M.Athans, R.P.Wishner and A.Bertoni, Suboptimal state estimation for continuous time non1inear system from discrete noisy measurements, IEEE Trans., AC-13, 504 (1968).
d ~----~r-----~r-----~b
f~------~------~~------~c
a
~------~------~~----~d
Fig.2. Converted diffusion flux on boundary line
b
Fig.l. Tetrahedron element
d'
e e a
bL------~
bL-----~
c
Fig.3. Converted diffusion flux on boundary surface
r -_ _ _ _ _ _~------~------~12
*
*
*
*
*
*
*
*
*
*
*
*
~----~~----~~------~8
Fig.4. Two-dimensional problem
*
*
• Filtered nodes
?: Filtered elements
*
*
0 Unfil tered nodes * Unfiltered elements
Fig.5. Application of filtering
652
Y. Oshima e t a l.
2 (1.2.4.5) (2.3.4.7) (2.5.6.7) (4.5 . 7.8 ) (2.4.5. 7 )
2
0.2.3. 6) 0 . 5. 6. 8 ) (3.6.7.8) 0.3.4.8 ) 0 . 3.6.8)
Fig.7. Divisio n of hexahe dron
Fig.6. Divisio n of Himeji city area
(a) Predic ted distrib ution at 5 oclock A.M.
(c) Predic ted distrib ution at 6 oclock A.M.
(b) Modifi ed distrib ution at 5 oclock A.M.
(d) Modifi ed distrib ution at 6 oclock A.M.
Fig.8. Predic ted and modifie d distrib utions of polluta nt concen tration
Prediction of pollutant concentration
Initial Conditions: S02 concentration 12 ppb Covariance of concentration 1.0 Diffusion coefficients Kx 200 m2/sec Ky = 200 m2/sec Kz = 280 m2 /sec Covariance of diffusion coefficients 10.0
ppb • Observed data
30
653
gComputational results
20
10
0 + - - - - - _ T - -_ _--~--~_4~_T--~--~--~_4--_T1--~--+_--r_~--_+------+_--~~
o
2
4
6
8
10
12
14
16
Fig.9. Observed concentration data and predic~ed concentrations at nodes surrounding monitoring station
18
20
h