NASA Contractor Report 178405 

ICASE REPORT NO. 87-72 

ICASE 

A METHOD OF BOUNDARY PARAMETER ESTIMATION FOR A 

TWO-DIMENSIONAL DIFFUSION SYSTEM UNDER NOISY OBSERVATIONS 

- - __ 

(NASA-CB-'l78405) A NEW HETHOD OF BOUNDABY 

PhRAtlETER ESTIMATION FOR A TWO-DZHENSZONAL 

DIPPUSIOB SYSTER UNDER HOISY OBSEBVAT ZOlS CSCL 128 Unclas 

~3/65 0111875 Final Report (NASA) 29 p 

P. Kojima 

Contract No. NAS1-18107 

November 1987 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, tlampton. Virgin4 23665 

Operated by the Universities Space Research Aseociatiop 
A KRTHOD OF BOUNDARY PA-R ESTIMATION FOR A TWO-DIMENSIONAL 

DIFFUSION SYSTEN UNDER NOISY OBSERVATIONS 

Y. Sunahara 

Faculty of Engineering and Design 

Kyoto Institute o€ Technology 

Kyoto, JAPAN 

F. Kojima 

Institute €or Computer Applications in Science and Engineering 

NASA Langley Research Center 

Hampton, VA 23665 

ABSTRACT 

The purpose of this paper is to establish a method for identifying unknown parameters involved in the boundary state of a class of diffusion systems under noisy observations. A mathematical model of the system dynamics is 

yiven by a two-dimensional diffusion equation, whose boundary condition is 

partly unknown diie to the existence of an unknown parameter. Noisy ohservations are made hv sensors allocated on the system boundary. Starting with the 

mathematical model mentioned above, an on-line parameter estimation algorithm 

is proposed wiLhin the framework of the maximum likelihood estimation. Existence OF the optlmal solution and related necessary conditions are discussed. Ry solving a local variatLon of the cost functional with respect to 

the perturbation of parameters, the estimation mechanism is proposed in a form 

of recursive computations. Finally, the €easibility of the estimator proposed 

here is demonstrated through results of d€gital simulation experiments. 

. 

Research was supported under the National Aeronautics and Space Administration 

under NASA Contract No. NASl-18107 while the second author was in residence at 

ICASE, NASA Langley Research Center, Hampton, VA 23665. 

i 
1. INTRODUCTION 

Recently, there has been much interest in the parameter identification 

problem for distributed parameter systems. For parameter identification problems for a class of distrthuted parameter systems, excellent survey papers 

have been piihlished by many researchers (e.g., Kuhrusly, 1976). Among them, 

the method €or parameter estimation incl.uding discussions on convergence 

properties o€ computer alqorithm has been developed by Banks and Lamm (19851, 

Yravaris and Seinfeld (1985), etc. Recently, feasible methods for estimat€ng 

unknown coefficients which depend on the state and spatial variables were repotted by the nlithors. In this paper, motivated by parameter identification 

problems of such parameters as heat conduction in chemical reactor or oil 

reservoir problems (e.<. , Chavent, 1978), we deal wIth the identif€cat€on of 

boundary parameters for a two-dimens€onal diffusion system under noisy ohservations. 

2. PROBLEM CONSIDERED 

[.et T and G be the time interval [O,tf] and system domain given by a 

bounded set in R , and let aC be the boundary of domain C. The boundary aG is characterized by C2-pararneterized curve decomposed into two 

parts such that 

2 

aG = aG"u aGK. (2.1) 

The system state u(t,x) at time t and at the location x = (xl,x*) is 

governed by heat diffusion equation: 

I 
-2- 

in Tx G 

toaetlier wtth the tnitLat and boundary conditions 

(2.2'3) 

(2.2b) 

(2.2c) 

(2.2d) 

22 22 where A = a /ax, + a /ax,, a/an denotes the outer normal. derfvative and 

6(6) is unknown parameter with a form of 

and where b(S! denotes the heat transfer coefficient. Since the heat 

transfer coeEficLent b has its physical value of positive definite, then 

0 F e(6) - < 1. It is apparent that, Cn the case where O(c) = 0, the boundary condition becomes the Dirichlet type while, in the case where e(€) = 1, 

we have the boundary condition OF Neumann type. 

N 

w 

Observation data are acquired through the sensors located on the boundary 

part aGK, as shown in Fig. 1. The observation mechanism is modeled by 
-3- 

‘/ 

x2 1 

0:Allocation of Sensors 

xl 0 

Fig. 1. Location of sensors on the boundary 

where y(t) = ly,(t),**-,y,!t)]’ and v(t) denotes the observation noise 

term modeled by an m-dimenstonal Srownian motion process. 

denotes the signal process defined by 

In Eq. (2.4), HI*] 

,(i) 

K where 

denotes the coordinates of sensor locations (zl , 

hi (i = 1,***, m) are assumed to be some constant, and where 

($1 $1). 
-4- 

The problem considered here is to flnd a method for identifying 

parameter H (t defined on a cU from information of the a priori 

hounct<jry known iti?lit g(S) and the noisy observation data {y(t))t>o. 

- 

3. THEORETICAL ReSULTS FOR PARAMETER ESTIMATION 

and p2 be meastires induced on the proc:esses y(t) and 1 Let u 

..# 

y(t) = v(t), respectively. Then, the Radon-Nikodym derivative is 

where 0 denotes the true function of unknown parameter. Associated with 

(3.1), the likelihood ratto functlonal is given by 

0 

The maximum likelihood estimate for the parameter 8 is the solution of 

where (3 denotes the admissi5le parameter class and 
-5- 

(See Raladrishnan, 1975, €or more details.) 9y using the gradient method (see 

Polak, 197 l), the opt Lmal solution can be obtained by the following recursive 

computations: 

(3 . Sa) 

where V, denotes the gradlent of the cost It(e). In order to compute 

V I (e(i)), 

each 0 : 

we have to SOIVP the following partial differential equation for 

Ot 

(€1 

-2- au(t') - Au(t,x) = 0 in T x G (3.6a) at 

with 

u(0,x) = u0(x;B (i)> on G (3.6b) 

This implies than in order to solve (3.3) massive computations are required. 

A possible method for avoiding this difficulty is to replace H[u(t;B)J in 

(3.2) by the statlonary val.iie H[G(O)] where u(x;e) = lim u(t,x;B) and 

where u(x;e) is a solution of 

- 

t+m - 

Aii(x) = 0 in G (3.7a) 
I 

-6- 

with houndary conditions 

on aGK. (3.7c) 

Thus, we introduce the Following cost functional: 

The Eollowiqg result states the relation between the likelihood ratio functional and the cost functional in this paper. 

Lemma 1: (See Sunahara, et al., 1984) The cost functional (3.8) satisfies the -- followiqg asyrnpJztt.: 

* 

Thus, our problem is to Find the optimal solution 8, satisfying 

(3.10) 

In this paper, we shall define the admlssible class of paraneters as follows: 
-7- 

[Definition] Le @? be the set of a1.l admissfhle parameters. Then, @J is 

saiA to be the admissLble parameter, i.e., 0 E (@, if 8 satisEies the 

€01 I ?wing proper t i es : 

The next results give the exfstence property for the optimal solution of 

(3. IO). 

Theorem 1: Suppose that 

Then, there exlsts at least one solution of (3.10) Eot A fixed t T with 

probability one. 

-- - 

--- 

The proof of Theorem 1 will be shown in Appendix I. 

The next results show the necessary cond€tion for the optimality of this 

paranetcr estimation problem. 

* 

Theorem 2: - tet 8 be the optimal solution of (3.10). Then the 

necessary condition for optimality is characterized by the following varia- 

- tionat inequalfty: 

-- 
-8- 

(3.11) 

* 

for Ye, Ot E 0 

-* * where uC0,) ts t'le solution of (3.7) correspondtng to 0 = et - and 

where ~(8~) Is the solritton of - the following Laplace equation with the 

rionhornogeiieoiis boundary condition, --- 

- 

* 

(3.12a) 

(3.12b) 

(3.12~) 

The proof of Theorem 2 will he shown in Appendix 2. 

By usfng the projection gradient method with the fixed step stze A, the 

computational algorlthn €or solving (3.10) can he conside,red as follows: 

(3.13a) 

(3.13b) 
-9- 

wlie re denotes the projection operator on @. Then, applying Theorem 1, 

(3. 13c) 

In the sequel, we propose the on-line parameter estimator. Let to and tf 

be the initial and the terminal time for parameter estimation. We choose the 

time stpp IC as 

tf - k I - 

n (3.14) 

and the tfmc interval (tg,tfl is divided into 

ti = to + ik for i = 1,2,*..,n-l. (3.15) 

I.le compute the reclirsive algorithm (3.13) at each time by using the 

observed data y(t;")). Considering the ffxed step size X in (3.13) as 

the time step k, the estimated sequence 8 can be obtained at each 

time 

tin) 

(*) for i = 1,2,.*.,n. ti 

[G(f)ln is characterized by the i= 1 Theorem 3: The estimated sequence 

(3.168) 
-10- 

for Y $ E 0 

i = 1,2,***,n-1 

- 

where u and p are the sol.utlons of 

(3.17) 

(3.16b) 
-1 1- 

2 

€or Y E Y (GI. JI 

4. ESTIMATION ALGORITHM 

4.1. Roundary Element - Approach --- 

In order to implement the proposed algorithm in the previous section, we 

adopt the boundary elenent method (REM) (see Rrehbia, 1978). Let v(x,xi) be 

the Fiindamrntnl solution satisfying 

where 6(x - xi) is the Dirac delta function. It is well-known that 

v(x,x') is given by 

where 

Ry using Green's Eormula and from (4.21, Eq. (3.7) can be rewitten by the 

follow€ng boundary integral equations: 

r(x,xi) denotes the distance between x = (x1,x2) and xi = (xl,x2). ii 
-1 2- 

where Ci denotes the constant such that 

i 

i 

1 for x E G 

l/? For x E aG 

0 for xi E F . 

i 2n c p(x ;e> i 
-1 3- 

In order to perform BEM, the boundary 36 is approximated by aE, which 

is decomposed into (n + m) elements, i.e., as illustrated in Fig. 2 

*2 

8% = y a'ZS, (4.6a) 

x1 0 

Fig. 2 Boundary elements and nodes of the system. 

(4 . 6c) 

In Fig. 2, (Pu (j))n and {PK (j) )j=l m denote the node of boundary elements 

EaGu (j) )j=l n 

(j) m 

"K 'j=l 

j=1 

and {8GK (j) Ijml m 

coincide with the sensor allocation as shown in Fig. 1. For 

respectively. It is noted that the nodes 
-1 4- 

For economy of notations, lrittoduce thr Following vector forms: 

Then, from (4.4) and (4.5), the botindary steady state is com~.~,A~-ri by the Following vector-matrix equations: 

(4.9) 

where A(0) and C are (n + m) x (n + m) matrices. 

4.2 Estimation Algorithm --- -- 

The numer€cal. procedure for the estimation algorithm in the previous section can be presented as follows: 

Step 0: Select an Initial parameter 8 at time t = to > 0 such that 

and set i = 0. 
-1 5- 

Step 1: Compute {<(Pu;G(i))} and {p(Pu;B by solving (4.8) and 

(4.9). 

Step 2: Compute the estimated parameter 8 at time t = to + (i + l)k 

by 

p+n) (1)) p) (j)) 

(*U (pU 

and 

(4.11a) 

(4.11b) 

Step 3: Sett€ng i by i + 1, go back to Step 1. 

5. NUMERICAL EXPERIMENTS 

Throughout numerical experiments, the system domain is given by a rectangle as illustrated in Fi2. 3. Figure 3 shows also that the boundary is 

divided into 20 elements, while ac, and aGK are decomposed into 5 and 

15 elements, i.e., n = 5 and m = 15, respectively. The boundary input g 

is set as 

R(5 ) = 500 on 3GK. (5.1) 
-16- 

Observations arc inane in the Form of 

(5.2a) 

= 0 €or j = 1,***,15, i = 0,1,2,*.*, (0) 

'i 

I 

I 

5. 0' - 

System Domain 

0 5.0 10. 0 

Fiq. 3. System domain G and its boundary aG 

where 

h = 3.1 for j = 1,***,15, j (5.2b) 
-1 7- 

arid wlwrc AV denotcbs ai1 Increment of tht- standard Brownian motion 

process qvnt>rntc4 by the Cnussf.?n distrihution. Tn Eq. (5.3), the system 

stnt~ (I( t,x;cr) is computed by adopting REM. 

(i) 

i 

Example 1: The unknown boundary parameter is preassigned as 8 = 0.2. By 0 

suhstltiiting known boundary input g and observation data {y (3) } 

(i = 1,2,***) into the proposed algorithm, starting with an initial 

parameter rit time to = 200 x k (k = l.O), the estimated parameter 

could be computed. Results of this example are demon- G(i)(s 1) i=O, 1,m.m 

strated in Fi%. 4. 

0 P 

- Iteration 0 

/ Iteration 600 

/ 

,,- 

Iteration 1200 

Iteration 1800 

Iteration 2500 

------do([): True Value 

t 

Fig. 4. Tnc ".;tirnsred parameter in Example 1. 

Example 2: In Example 2, we deal with the case where the unknown parameter is 

a space-variable function. The unknown parameter and simulation results are 

illustrated In Fig. 5. 
it(€). L 

1.0 

0.5 

0 

-do(f): True Value 

-18- 

-Iteration 0 

-Iteration 500 

-Iteration 1000 

- Iteration 3000 

- Iteration 6000 

P 

;ig. 5. The estimared pa:aneter in Example 2. 

6 0 CONCLUSIONS 

A method €or identifying the boundary condition has been developed in 

this paper For a distrihiited system whose dynamics are governed by a twodimensional heat diffusion equation. Rased on the concept of maximum likelihood estimate, the problem was converted into the optimal control problem. A 

difficulty arises in the computational burden involved in computing the gradient of the associated likelihood ratio functional. In this paper, to avoid 

this difficulty, the steady state model was introduced. The validity of the 

estinator proposed here was demonstrated through numerical experiments. 
-19- 

REFERENCES 

Rala,lrishnan, A. V. (1975). "Identification and stochastic control of a class 

of distributed parameter systems with boundary noise," Lecture Notes in Economics arid ?fathematical Systems, Vol. 107, Springer-Verlag. -------------- 

Banks, 11. T. and P. 0. Lamm (1985). "Estimation of variable coefficients in 

parabolic distributed system," TEEE Trans --., AC Vol. AC-3CI, pp. 386-398. 

Rrebbia, C. A. f1978). The Boundary Element Method €or Engineers, Pentech 

Press, London. 

Chavent, C. (1978). "About the identification and modeling of miscible or 

immiscible dtsplacements in porous media," Lecture Notes in Control and 

Tnformat ion Sciences, Vol. 1, Springer-Verlag, New York. 

Ko.jirna, F. I Y. Srinahara and <:h. AI hara (1985). "Parameter estimation for distributer? systrarns with obstacle under noisy observatton," Proc. 7th TFAC 

-I---- 

Syinposiim - - --- on l____l-_-- Tdentification and System Parameter Estimation. - 

Kravsris, $. and .J. H. Seinfeld (1985). "Identification of parameters in distributed parameter systems hy regularization," SIAM J. Control Optim., Vol. 

23, pp. 217-241. 

Kubrusly, C. S. (1977). "Distributed parameter system identification. A 

s:irvey," --- Int. .J. Control, - 1'01. 25, pp. 509-535. 
-? 0- 

Pol.ak, E. (1971). Computational Methds ---- in Optimization, Academic Press. 

Sunahar?, Y., Sh. Aihara, and F. Kojima (1984). "A method for spatial domain 

identification o€ distrlhuted parameter systems under noisy observations ," 

Proc. 9th IFAC World Conpress. 

Siinahara, Y., Sh. Aihara, and F. Kojima (1986). "A method €or parameter estimation of a class of nonlinear dtstributed parameter systems under noisy observations," Automatica, Vol. 22, pp. 727-732. 

, 
-2 1- 

Appendix 1: The Proof of Theorem 1 

, 

1,ct (0 ] €(.@ he a minimizing sequence such that 

n 

Noting that 

and 

Using the compactness nethod, it can be shown that, under the condition (C-I), 

I 

u(8n) + ~(8) strongly in H*(G). (A.6) 

Nence, we have 
-22- 

E{ IJt(en) - Jt(0)JI + 0. (A.7) 

- * 

The proof has been completed. t' 

From (A.l) and (A.5), we may Find 0 = 0 

, 
-23- 

Appendix 2. The Proof of Theorem 2 

if : k ;;>t imjl solution exists, the followlng variational inequaldenotes the gradient of Jt(0). In order to obtain V,.Jt(9), ‘e where 

let tis derive the local variation of cost functional (3.8). Let OA he the 

!wrturhed parameter of 8 such that 

where h is a positive constant and de([) is a glven real valued function. Set RS 

Applying (3.8) to the mean valtie theorem, (A.lO) becomes 

6Jt(8) = -H[6U(8)]’ {, y( t) - H[Y(O)]). 

\ 

On the other hand, 6u(x;B) ts a solution of 

(A.ll) 

(A. 12) 

sT;<x) = o in G (Ao13a) 
-24- 

I Introdlice the so-called “ad joint system” of which state 

tion o€ 

Ap(x) = 0 in G 

with 

(A. 13b) 

. 

(A. 13c) 

is a solu- 

(A. 14a) 

(A. 14b) 

(A. 144 

where Irk[*] denotes the adjotnt operator of H[*l. By virtue of Green’s 

formula, from (4.13) and (A.14), we have 

(A. 15) 

From (A.10) and (A.13), 6Jt(e) is computed by 
-2 5- 

. Thus, we obtain 

1 VeJt(8)-S0 = lim 7 6Jt(8) 

A+ n 

(A. 16) 

(A. 17) 

* * 

t By setting 8 = 0 and 

tional inequality (3.11). 

63 = 8 - et in (A.171, we can obtain the varia- 

\ 
Report Documentation Page 

1. Report No. 2. Government Accession No 

___~__- ~- ~ ~___ 

7. Authorls) 

Y. Sunahara and F. Kojima 

UnclassiEied 

__-_ __ -- --. .- --___ 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 

Mail Stop 132C, NASA Langley Research Center 

Hampton, VA 23665-5225 

and Engineering 

Unclassified 27 A0 3 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Langley Research Center 

Hamp ton, VA 2366 5-5225 

15. Supplementary Notes 

Langley Technical Monitor: 

Richard W. Barnwell 

Final Report 

3. Recipient's Catalog No. 

5. Report Date 

Novembe r 1987 

6. Performing Organization Code 

8. Performing Organization Report No. 

87-72 

10. Work Unit No. 

505-90-21-01 

11. Contract or Grant No. 

NAS1-18107 

13. Type of Report and Period Covered 

Contractor Report 

14. Sponsoring Agency Code 

Submitted to the Proc. of the 10th 

IFAC World Congress on Automatic 

Control, Munich, FRG, 1987 

16. Abstract 

The purpose of this paper is to establish a method for identifying unknown 

parameters involved in the boundary state of a class of diffusion systems under 

noisy observations. A mathematical model of the system dynamics is given by a 

two-dimensional diffusion equation, whose boundary condition is partly unknown due 

to the existence o€ an unknown parameter. Noisy observations are made by sensors 

a1 located on the system boundary. Starting with the mathematical model mentioned 

above, an on-line parameter estimation algorithm is proposed within the Eramework 

of the maximum likelihood estimation. Existence of the optimal solution and related necessary conditions are discussed. By solving a local variation of the 

cost functional with respect to the perturbation of parameters, the estimation 

mechanism is proposed in a ferm of recursive computations. Finally, the feasibility of the estimator proposed here is 

simulation experiments. 

__-___ 17. Key Words (Suggested by Author(s)) 

identification, parameter estimation, 

distributed parameter systems, 

stochastic systems 

demonstrated through results of digital 

____ 

18. Distribution Statement 

65 - Statistics and Probability 

66 - Systems Analysis 

Un class i f ied - un 1 imi ted 

=--No. of pages 122. Price