- Email: [email protected]
A Method for Spatial Domain Identification of Distributed Parameter Systems Under Noisy Observations
Cnp\right © IF.\( : ~'th I lil"lIlll.d " ·"Ilt! Budape .. t . 111I1l.l{. ln . 1~ ' :-; -4
("Il ~Il""
A METHOD FOR SPATIAL DOMAIN IDENTIFICATION OF DISTRIBUTED PARAMETER SYSTEMS UNDER NOISY OBSERVATIONS Y. Sunahara, Sh. Aihara and F. Kojima
Abstract : Motivated by oil reservoir problems, this paper is concerned with the theoretical and computation aspects of a method for identifying the spatial domain of distributed parameter systems under noisy observations . The system model is given by a partial differential equation of parabolic type derived by a known distributed input . A mathematical model with an additive noise term of the observation mechanism is also given . Based on the concept of the maximum likelihood estimate, the estimation algorithm is presented in a form of the recursive computation . The main body of theoretical aspects in this parer is convergency properties of the estimation algorithm. Both the feasibi ity and the validity of the method presented here are discussed including results of digital simulation experiments . K7ywords : Ident~fication ; par~meter .est~mation; partial differential equatlons ; stochastlc systems; maXlmum llkellhood estlmates . 1.
INTRODUCTION AND STATEMENT OF THE PROBLEM
u(o.x) = uO(x ) on Go (1b) u(t.xo (O) = O. on ]O . T[ . xo(~) eaG o and
There has recently been much practical interest shown in lroblems of identifying the spatial domain 0 a class of distributed parameter systems or finding the optimal shape of structural systems . For example, motivated by oil reservoir problems, identification and mode ling problems of miscible or immiscible displacement in porous media have been studied , where an identification problem has been formulated as a control problem (Chavent , 1978 ). On the other hand , in optimal shape design problems , several excellent studles have been reported . (Pironneau, 1977 , 1982 and references cited in those papers)
~e[0.1J
where x=(x 1 .x 2 ) and naturally 6= a2~ axf+ a~ ax~. In (la). the distributed input f(x) is known and this means an added pressure on the known s ub-region 0 1 (c Go ) in the unknown region Go for the purpose of estimating the region Go ' Reactions of the system are observed with the background noise v(t) . Thus . setting m-dimensional observation mechanisms . observations are given by
The principal line of theoretical attack in the problem mentioned above is to convert the problem into the optimal control where the control variable represents an unknown spatial domain (Begis, et. aI., 1975 and Delfour, et al , 1983) Feasible algorithms for finding the optimal shape were also proposed by several authors where a class of partial differential equations of elliptic type was mainly considered as a system model of control problems
y(t) = rtH[u(s:G )]ds JO 0
f(x)
+
v(t )
( 2)
where y(t)=[Y1(t)'Y2(t) . · •• Ym(t)]' and v(t) is supposed to be an m-dimensional standard Brownian motion process . In (2). for arbitrary functions ~. the operator H is defined by
H[~]
In this paper , considering the fact that, in exploring problems of underground resources, the sytem state u (t , x) expresses frequently the pressure or the density of the media under consideration , a very simplified model of the system is given by the diffusion equation with its unknown time invariant domain Goc R(2). To fix the idea. we shall first describe the mathematical model of the system considered . au(t.x) _ 6U(t.X) at
(lc)
=
[~w
hI ( x )~( x)dx . 1
~ Wh2 ( x )~( x ) dx . 2
~ W hm ( x)~(x ) dx] '
(3)
m
and
The problem is to find the region Go provided that the input f(x) and the observed data y(t) are obtained respectively in DIand W. In Section 2. the problem at hand is concretely stated in a form of estimating the unknown xo(~)-curve so as to minimize the likelihood ratio functional. A method is presented in
in]O.T[XG o ' (la)
with the initial and boundary conditions
649
Y. Sunahara, Sh. Aihara and F. Koj ima
650
Section 3, for estimating xo ( ~) based on noisy observation data, In Section 4. convergence properties of the estimation procedure are discussed by introducing the concept of stochastic stability, Section 5 is devoted to showing results of digital simulation experiments, For convenience of present descriptions. the function spaces used here are listed below : L2 (G) : a space of square integrable functions on GcR(2). C= (G): a space of infinitely differentiable functions on G . (G): a space of functions such that { cp I
cO'
cp e C= (G). cp with compact support in G } C2 (G): a space_of twic~ differentiable functions on G where G is the closure of G H6(G): a space of functions such that { cp I
v
sets. the maximum likelihood estimate with respect to the domain G is a solution to
Since massive computations are required to solve (7) with (6). we replace the system state u(t.x ) by its stationary state u ( x ) =lim _ =u(t. x ). where x E G, Thus. we write t for (5)
X [t.u(G)]
exp[H[u (G)]'y( t ) - ~ tH[u (G)] ' H[u(G )]] (8)
where u satisfies - t.u(x)
f(x)
in G •
( 9a )
u(x(~) )
0
on aG ,
( 9b )
The problem is to find G* in such a way that for t> 0 . (10) where, instead of ( 6).
2,
MAXIMUM LIKELIHOOD ESTIMATE
The following three assumptions are made : (A-I) there exists two open bounded regions Dl and D2 in R(2 ) such that Dl cGocD Z : (A-2) feC Q (D 1 ) : (A-3) u eL 2(D ) , O 2 Re~arding the observation mechanism given by (2).the following conditions are set : (C- 1) Ui::omWic Dl and WinWj= ~ (i*j. i. j=1. 2.' • • . m) : (C-2) hieCQ(W i ) ( i=1.2,··· .m) Let (a.F.p) be a basic probability space and let ~1 and ~Z be the measures induced by the processes y(t) and y(t)=v(t).respectively, Then, the Radon-Nikodym derivative is given by
~ (t.u (t :G )) = exp[ ~~ U/"2
itH[u (t:G )]'dy(s)
0
0
0
- -,l,..-J\[u ( s:G )], H[u ( s:G ) ]ds]. ( 4) .. 'JO 0 0 We define the likelihood functional by A[t.u(t:G )]
The following theorem states an asymptotic property of Jt(G ), [Theorem-I] The cost functional Jt (G) converges to Jt (G) in the following sense lim E{ I Jt (G) - Jt(G ) t-=
(12)
[Theorem-2] With the conditions (C- 1 ) to (C-4) , there exists at least one solution to (10) for a fixed we a ( See Chenais 1975 ) 3,
ESTIMATION ALGORITHM
In this section. we shall develop a feasible algorithm for finding a solution to ( 10 ), To do this. let GA be the domain such that aGA = {x (~) + AOt ( ~ ) n ( ~ ) I ~ e [0. 1] } . (13) where n(~ ) ( nl ( ~) , n 2 ( ~ )) and Ot ( ~ ) are respectively a normal vector and a given real valued function at x ( ~ ) = (x l ( ~ ).x 2 ( ~ ))' Along the result by Pironneau ( 1977 . 1982 ). the following lemma can be easily obtained by using (8) . (11) and ( 13 ), [Lemma] Define
( 5)
where Gc R( 2) and introduce
f} = 0
OJt(G)~Jt(Gl) - Jt (G) Then.
O\(G) = A
r
JaG
Ot (~ )oii(~ ( ~ )) op(x(n )dl an
on
+ o(A).
Let G be a class of admissible sets, The domain G is said to be admissible . i. e,. G e G if the following conditions hold: (C-3 ) G is a bounded open set in R(2 ) such that D1 cGcD 2 : (C-4) G is the interior region of a C2-parameterized curve given by {x( ~ ) I ~ e [0. 1] } , With the definition of a class of admissible
where p(x ) is the solution to (vp ( x )) ' vcp(x )dx = -H[cp] , {~ G
~
(15)
i
- H[ii(G ) ] } .
( 16 )
for any cp e H6 (G). dl~ denotes the line element and ll..:.2 ~ on
v(' )
(18)
Spatial Doma in Identifi cation of Distribut ed Parame t e r Sys t ems
We shall implement an algorithm for finding - (G* )=0. As stated the domain G* such that 8~ in Lemma . 8~ (G) is obtained only through the boundary curve x ( e ) of the domain G. Then. for l>O. set as _ aJ (G ) xn+1 (e )=xn (e)-l 'Qa n n (e ) ( 19a ) ( 19b ) where Gn is the domain obtained by the n-th recursive computation and the boundary of Gn is denoted by aGn={xn (e ) leE [0 . 1] }. In (19b ). the operator PG is a projection on G. Denoting the boundary integral on aG by <• . • >oG . i . e . . for any ~ and ~.
we may write 8J .(G.OI ) =
aj .( G)
-~p ( x ) = -h(x ) ' { ~(t)-H[u ( Gt ) ] } for x E Gt with the boundary condition for
e E [ 0. 1]
( 25a) ( 25b )
where u(Gt ) is the solution to ( 24 ) and h( x ) =[h (x ) . h (x ) . . • • . hm (x )] ·. By solving 1 2 ( 24 ) and ( 25 ) . a pair of the solutions {u (Gt ). p(G t ) } is obtained . We (Step-3 ) Es timat ed domain co ns true tion shall set the perturbation parameter l in ( 23 ) as l=k . Then. the perturbation formula ( 23 ) is rewritten as
-
-
xt+k (e)=xt (e )-k
aut ~t (e ) ) an
( 26a) ( 26b)
< 01
•
aj .( G) --aa-
( 21)
>aG
Thus . comparing ( 21 ) with ( 15 ) . we obtain
--aa-
651
a ~ ( x ( e ))
on
op( x(e )) on
(22)
Hence (19 ) is _ o~ ( xn ( e )) op ( xn (e )) xn +1 (e )=x n (e)-l an an n(e ) ( 23a) ( 23b ) Based on ( 23 ). the following procedure to perform_computer implementation is stated : Let Gt be the estimated domain of Goat time t. We denote the boundary of Gt by aGt= { ~t ( e ) I eE[O.l] }. (Step-I ) Preassignment of an initially estimated domain At time to we select an initially estimated curve ~t ( e ) E G such that
o
~t ( O )· ~t ( 1 ) and construct the estimated doo _ 0 main Gt with the boundary oG t = { ~t ( e ) leE
o
0
0
[0 . 1]}. Let t be the discrete-time index t=tO+ik . i=O . 1. 2.3 . • . • • • . where k is the partitioned time . By using a newly observed data at each time t . recursive calculations for computing the estimated boundary curve associated with the domain Gt ( t = to . to+k. t o+2k . • • • ) are performed along Step-2 and Step-3.
The estimated domain Gt +k is determined by the curve oGt+k= { ~t+k ( e ) leE [0 . 1] } . Convergence properties of Gt are examined along the theoritical aspect which will be stated in the following section . Remark: The projector PG in ( 26b ) implies that the boundary curve ~t+k ( e ) is given by joining together the curve ;t+k (e ) and the boundary curve of D . Thus . it may reasonably be as2 sumed that the boundary curve ~t ( e ) is modified so that ;t (e ) may become C2-curve . 4.
CONVERGENCE PROPERTIES OF ESTIMATION ALGORITHM
It is difficult to examine asymptotic behaviors of the estimated domain Gt because the domain estimate i s performed by the estimated boundary curve { ~ ( e ) leE [0 . 1] } based on the solutions to ( 24 ) and ( 25 ). A possible method for investigating asymptotic behaviors is an evaluation of the cost functional jt (Gt ). Starting with ( 8 ) and ( 11 ) and noting that . from (14 ) . jt+k (Gt+k )-jt (Gt +k ) ( 29 ) 8jt (Gt )
jt+k (Gt+k )-jt (Gt ) +
we obtain
(Step - 2 ) Steady s tate solution of th e sy s t em and its adjoint one . Using the estimated domain Gt for G in ( 9 ). we write -.:.iI( x ) = f (x ) for x E Gt (t=tO+ik ) ( 24a ) with the boundary condition for
e E [ 0.1]
( 24b )
On the other hand . the adjoint equation ( 16 ) is expressed by
From Lemma and ( 20 ) . it follows that 8J (G ) = -k(- ( O~ ( x ( e )) ap(x(f ))) 2dl t t JaG on an t
+
o(k) .
e
( 31 )
Y. Sunahara, Sh . Aihara and F. Kojima
652
L W(t . q )
Then. from ( 30 ) and ( 31).mu1tiplying both sides of ( 30) by 1/ k. integrating the result over [to ' T] and letting k-+ O. we have -8 (t )dt + ~ H[u(G t )] '{ ~ y(t) - H[u(t . Go )] }dt - ~ H[u(Gt)] ' dv(t ) (32 )
dJ(t)
with
xexp{
r_
J,oG
( ou(~(e» on
op(x(f) ) )2dl
on
From (P-l) to (P-3 ). it follows that
(34 )
e
t
The following proposition states ~n asymptotic behavior of the solution pr~cess J (t ) of ( 32 ): [Proposition] The solution J (t ) to ( 32 ) satisfies the following equation. i,e , . for any ~>O . lim P[ 1 J (t ) t-+ 00
Proof
+
~ y(t) ' y(t ) 1 2> d 2t
Set as
q(t) = J (t ) + ~ y(t) ' y(t ) , 2t
= 0, (35) (36 )
<0
-28 (t)q (t )
if q(t).;, 0
-8(t )dt+y(t )dt+ry(t )' dv (t)
L
( 42 )
W(t .q ) < 0 .
( 44 )
See Has ' minskii
Consequently. for any ~ > O . 1980 . for more details ) lim P[ t-+ 00
1
q (t )
2 >~ ] = 0 .
1
( 45 )
Thus . the proof has been completed . ~onvergence property of the estimated domain Gt to the true domain Go is shown by the following theorem .
~ > O,
[Theorem-3] For any holds :
Using Lemma. the q (t )-process satisfies dq(t )
(41 )
-I y( t ) 1 (q (t) )2+ 2y (t )q(t )- 1 y( t ) 1 = - 1 y(t ) 1 {( q(t ) )2_ ,2~ i ~ l , q (tH } :o; O. ( 43 ) Hence . if q (t)';' 0
where O(t)=
f,oo;. (s) 1 ds }. t
and
!J
{-28(t)q(t) + 2y (t)q (t )
:0;
- 1 y ( t ) 1 (q ( t) ) 2_ 1 y ( t ) 1 }
lim P[ t-+ 00
( 37a )
the following limit
-
11
H[u(Gt )] - H[u (Go )]
11
(m»d=O.
R
(46 )
with -
1
'
J (t O) + - 2 y(t O) y(t O) . ( 37b ) 2t
Proof :Noting that lim v~t ) = 0 t-+ 00
where y (t ) and ry (t ) are a scalar and an m-dimensional process respectively given by y(t) = -~ [~y ( t ) - H[U(Gt ) ]] ' [~ y(t ) and
- H[u(t;G )]] + ~2 o 2t
( 39 )
lim t-+oo
~
lim t-+ 00
q(t)
0 ;
(P-2)
if q(t)=O . then 8 (t )=0. y(t )= ~ 2t2 and ry(t)=O
-
<
00
= 1
]
and
to
P{
L ry (t) , ry (t )dt < 00
t-+
00
=
Associated with (P-3). define W(t.q ) = { (q(t) )2 + LOC>I y(s )
t
oo
1
ds
t
ry (s) ' ry(s )ds }exp{Loo, t
00
respectively . From ( 35 ). ( 49 ) and ( 50 ). the equality ( 46 ) holds . DIGITAL SIMULATION STUDIES
Exampl e -I : The distributed system input f (x )
to
+L
t ........
+ ~ H[u(Gt )]'H[u(G t )] } w, p. l (50)
5. J
00
-
lim J (t ) = lim { -H[u(Gt )]'H[u(G o )]
if q(t)';' O. then 8(t » 0
P[ LOC>I y ( t) 1 dt
( 49)
w. p. 1
and
(P-l)
(P-3)
j(tH[u (s;G )]ds 0 0
it is easy to show that
With (C-1) and (C-2 ) . it can be shown that the q (t)-process has the following properties : ~
( 47 )
and
( 38 )
ry (t ) = ~ { ~ y(t) - H[u(Gt )] }.
w. p.1
y (s ) 1
ds } . ( 40)
and write the differential generator by L associated with (37 ) , we obtain
in (1) is given by f(x ) {=-2 forxEDicD1 ( 51a ) = 0 for x E D1 and -2 :0; f (x ) :0; 0 for xED 1 - Di ' (51b ) Letting m=l . i . e .• Wl=D l in ( 2). observations are made in the form of dy (t ) =
(j[ h1 (x )u(t.x )dx)dt + dv (t ) D1
(52)
653
Spa tial DomaiR Identificati on of Distributed Pa rame t e r Systems
with y(O)=O . where for x e Di . h (x)f=h (53a) 1 Lo for x e Dl and Os; h 1 ( x ) s; h forxeD1-Di . ( 53b) (Preassignment of Go and D1 ) In Example-I . for the purpose of simulation experiments. the true spatial domain Go which is to be identified and the observation domain Dl were preassigned as shown in Fig . 1. The domain Go was decomposed by 101 triangles of subregions with 64 nodes . curve of true domain Go
X
la'
dye tl/ d t
100
50
o
50
100
Time
Fig . 2. Sample run of the scalar observation dy (t )/dt in Example-l 150
(52). respectively. obtain the solutions u and p to ( 24) and (25 ). Step-3 : By substituting the solutions and p obtained in Step-2 into ( 26 ). the k-th iterated curve of the estimated boundary ~t +k (E ) is computed . 0 St ep-4 Construct the estimated domain Gt +k
u
100
50
o
o
50
10 0
150
200
Xl
Fig . 1. Unknown domain Go ' given region Dl and their triangulations in Example-l (Observation data ) In simulation experiments. the observation data dy (t ) were generated by ( 52 ) and ( 53 ) . where the u( t . x )-process was obtained by applying the finite element method to ( 1 ) with the input given by ( 51 ) and the s tandard Brownian process v(t ) generated by a usual way. A sample run of the observation process is shown in Fig. 2. (Domain identificaion ) Based on the computer algorithm stated in the previous section . the numerical procedure is as follows : St ep-l Preas sign the initial curve { ~t (E)} as shown in Fig . 3. 0 St ep-2 With the input f and the observation data . i.e .. the output data. given by ( 51 ) and
based on the boundary curve derived in Step-3 . Letting k=I.2 . • • • . Steps ( 2 ) to ( 4) give a forwardly recurre~t algorithm to obtain the running estimate Gt with the input f ( x ) and the output dy(t ) as a set of the given data . Figure 3 shows the evolution of estimated boundary curves .
Example-2
The distributed input f (x l is
given by
f(x ) [ :
-40 0
for x e Di for x e Dl
and -40 S; f ( x) S; 0 for x e D1- Di ( 54b ) The observation mechanism is set with fiftyone observation domains. loll' 101 2 , •• 101 51 , (Preassignment of Go and Dl in Example-2 ) Regarding the true domain Go and the preassigned domain D1 . the same decompositions as in Example-I were adopted. These are illustrated in
1500
15 .0 iteration 1580 curve
,, ,,
lOO
,
··· ./
5.0
,,'/
,// ./
o Fig . 3.
( 54a )
15.0 5.0 10.0 20.0 Estimated boundary of unknown domain Go in Example-l
Y. Sunahara, Sh. Aihara and F. Kojima
654
Chavent . G ( 1978) . About the identification and modeling of miscible or immisicible displacements in porous media In Lec ture Not es in Control and Information Sci ences, Vol 1 . pp . 196-220 . Springer. New York . Chenais . D (1975 ) On the existence of a solution in a domain identification problem j. Math. Anal. Appl. 52. pp 189-219 Delfour . M.. G Payre and J-P Zolesio ( 1983 ) Optimal design of a minimum weight thermal diffuser with constraint on the output ther mal power flux . Appl. Math. Optim. 9 . pp , 225-262 St ochas ti c StabiliHas' minskii . R. Z. ( 1980) . ty of Differ ential Equations , Sijthoff and Noordhoff . Alphen aan den Rijn Pironneau . 0 (1977 ) Variational methods for the numerical solutions of free boundary problems and optimum design problems . In Control Th eory of Systems Governed by Partial Differ ential Equations, 209 . Academic Press . New York . Pironneau . 0 (1982 ). Optimal shape design for elliptic systems . In Lec ture Not es in Control and Information Sci ences , Vo1. 38 . pp , 42-66 . Springer New York . ACKNOWLEDGMENTS The authors wish to e xpress their thanks to Messrs . Teruo OKA and Toshiyuki MORIOKA for their helpful assistance .
Fig
4 . Each observation domain \'\ (i=1.2.· . 51 ) was given by a triangle as shown in Fig 5 Using the 51-dimensional observation mechanism . a sample run of the estimated domain is shown in Fig 5. iterated by the same procedure as described in Example-I . CONCLUSIONS A feasible method has been presented in this paper. for identifying the boundary curve of the unknown domain of distributed parameter systems under noisy observations The following three aspects in this paper should be emphasized as concluding remarks ( 1 ) Math ematical model Recognizing the fact that dynamic behaviors of many real physical systems are . in fact. distributed. a mathematical model is described by partial differential equations of parabolic type not elliptic . Furthermore. noisy observations are taken into account from viewpoints of practical applications . (2 ) Recursive algorithm for computing th e unknown spatial domain By using current information of the input-output data . the unknown spatial domain can be estimated by recursive computations . (3 ) Us e of th e stationary state information of th e system A difficulty arises in computing recursively the unknown spatial domain because of increase of the observation data . Consequently . use is made of the steady state solution of the system for the purpose of circumventing massive computations . Studies on convergency properties described in Section 4 support the feasibity and the validity of the proposed method . Turning our eyes to context of shape design by Pironneau (1977 . 1982 ) and Begis . et al ( 1975 ). it can be understood that the idea involved in this paper bridges between the method of the optimal shape design and that of maximum likelihood estimate
domain Go region Dl
50
REFERENCES Begis . D. and R Glowinski (1975 ). Application de la method des e lements fini a l ' approximation d ' un probleme de domain optimal . Methodes de resolution des problemes approches Appl. Math . Opti m. 2 . pp . 130 - 169 .
o
50
100
150
iteration 7500 15 .0
true curve
,
:' 10.0 :,
,,
5 .0
'-'
o
5.0
Fig. 5.
100
200
X,
Fig . 4. Unknown domain Go ' given region D1 and their triangulations in Example-2
15.0
200
Xl
Estimated boundary of unknown domain Go in Example-2
("Il ~Il""
A METHOD FOR SPATIAL DOMAIN IDENTIFICATION OF DISTRIBUTED PARAMETER SYSTEMS UNDER NOISY OBSERVATIONS Y. Sunahara, Sh. Aihara and F. Kojima
Abstract : Motivated by oil reservoir problems, this paper is concerned with the theoretical and computation aspects of a method for identifying the spatial domain of distributed parameter systems under noisy observations . The system model is given by a partial differential equation of parabolic type derived by a known distributed input . A mathematical model with an additive noise term of the observation mechanism is also given . Based on the concept of the maximum likelihood estimate, the estimation algorithm is presented in a form of the recursive computation . The main body of theoretical aspects in this parer is convergency properties of the estimation algorithm. Both the feasibi ity and the validity of the method presented here are discussed including results of digital simulation experiments . K7ywords : Ident~fication ; par~meter .est~mation; partial differential equatlons ; stochastlc systems; maXlmum llkellhood estlmates . 1.
INTRODUCTION AND STATEMENT OF THE PROBLEM
u(o.x) = uO(x ) on Go (1b) u(t.xo (O) = O. on ]O . T[ . xo(~) eaG o and
There has recently been much practical interest shown in lroblems of identifying the spatial domain 0 a class of distributed parameter systems or finding the optimal shape of structural systems . For example, motivated by oil reservoir problems, identification and mode ling problems of miscible or immiscible displacement in porous media have been studied , where an identification problem has been formulated as a control problem (Chavent , 1978 ). On the other hand , in optimal shape design problems , several excellent studles have been reported . (Pironneau, 1977 , 1982 and references cited in those papers)
~e[0.1J
where x=(x 1 .x 2 ) and naturally 6= a2~ axf+ a~ ax~. In (la). the distributed input f(x) is known and this means an added pressure on the known s ub-region 0 1 (c Go ) in the unknown region Go for the purpose of estimating the region Go ' Reactions of the system are observed with the background noise v(t) . Thus . setting m-dimensional observation mechanisms . observations are given by
The principal line of theoretical attack in the problem mentioned above is to convert the problem into the optimal control where the control variable represents an unknown spatial domain (Begis, et. aI., 1975 and Delfour, et al , 1983) Feasible algorithms for finding the optimal shape were also proposed by several authors where a class of partial differential equations of elliptic type was mainly considered as a system model of control problems
y(t) = rtH[u(s:G )]ds JO 0
f(x)
+
v(t )
( 2)
where y(t)=[Y1(t)'Y2(t) . · •• Ym(t)]' and v(t) is supposed to be an m-dimensional standard Brownian motion process . In (2). for arbitrary functions ~. the operator H is defined by
H[~]
In this paper , considering the fact that, in exploring problems of underground resources, the sytem state u (t , x) expresses frequently the pressure or the density of the media under consideration , a very simplified model of the system is given by the diffusion equation with its unknown time invariant domain Goc R(2). To fix the idea. we shall first describe the mathematical model of the system considered . au(t.x) _ 6U(t.X) at
(lc)
=
[~w
hI ( x )~( x)dx . 1
~ Wh2 ( x )~( x ) dx . 2
~ W hm ( x)~(x ) dx] '
(3)
m
and
The problem is to find the region Go provided that the input f(x) and the observed data y(t) are obtained respectively in DIand W. In Section 2. the problem at hand is concretely stated in a form of estimating the unknown xo(~)-curve so as to minimize the likelihood ratio functional. A method is presented in
in]O.T[XG o ' (la)
with the initial and boundary conditions
649
Y. Sunahara, Sh. Aihara and F. Koj ima
650
Section 3, for estimating xo ( ~) based on noisy observation data, In Section 4. convergence properties of the estimation procedure are discussed by introducing the concept of stochastic stability, Section 5 is devoted to showing results of digital simulation experiments, For convenience of present descriptions. the function spaces used here are listed below : L2 (G) : a space of square integrable functions on GcR(2). C= (G): a space of infinitely differentiable functions on G . (G): a space of functions such that { cp I
cO'
cp e C= (G). cp with compact support in G } C2 (G): a space_of twic~ differentiable functions on G where G is the closure of G H6(G): a space of functions such that { cp I
v
sets. the maximum likelihood estimate with respect to the domain G is a solution to
Since massive computations are required to solve (7) with (6). we replace the system state u(t.x ) by its stationary state u ( x ) =lim _ =u(t. x ). where x E G, Thus. we write t for (5)
X [t.u(G)]
exp[H[u (G)]'y( t ) - ~ tH[u (G)] ' H[u(G )]] (8)
where u satisfies - t.u(x)
f(x)
in G •
( 9a )
u(x(~) )
0
on aG ,
( 9b )
The problem is to find G* in such a way that for t> 0 . (10) where, instead of ( 6).
2,
MAXIMUM LIKELIHOOD ESTIMATE
The following three assumptions are made : (A-I) there exists two open bounded regions Dl and D2 in R(2 ) such that Dl cGocD Z : (A-2) feC Q (D 1 ) : (A-3) u eL 2(D ) , O 2 Re~arding the observation mechanism given by (2).the following conditions are set : (C- 1) Ui::omWic Dl and WinWj= ~ (i*j. i. j=1. 2.' • • . m) : (C-2) hieCQ(W i ) ( i=1.2,··· .m) Let (a.F.p) be a basic probability space and let ~1 and ~Z be the measures induced by the processes y(t) and y(t)=v(t).respectively, Then, the Radon-Nikodym derivative is given by
~ (t.u (t :G )) = exp[ ~~ U/"2
itH[u (t:G )]'dy(s)
0
0
0
- -,l,..-J\[u ( s:G )], H[u ( s:G ) ]ds]. ( 4) .. 'JO 0 0 We define the likelihood functional by A[t.u(t:G )]
The following theorem states an asymptotic property of Jt(G ), [Theorem-I] The cost functional Jt (G) converges to Jt (G) in the following sense lim E{ I Jt (G) - Jt(G ) t-=
(12)
[Theorem-2] With the conditions (C- 1 ) to (C-4) , there exists at least one solution to (10) for a fixed we a ( See Chenais 1975 ) 3,
ESTIMATION ALGORITHM
In this section. we shall develop a feasible algorithm for finding a solution to ( 10 ), To do this. let GA be the domain such that aGA = {x (~) + AOt ( ~ ) n ( ~ ) I ~ e [0. 1] } . (13) where n(~ ) ( nl ( ~) , n 2 ( ~ )) and Ot ( ~ ) are respectively a normal vector and a given real valued function at x ( ~ ) = (x l ( ~ ).x 2 ( ~ ))' Along the result by Pironneau ( 1977 . 1982 ). the following lemma can be easily obtained by using (8) . (11) and ( 13 ), [Lemma] Define
( 5)
where Gc R( 2) and introduce
f} = 0
OJt(G)~Jt(Gl) - Jt (G) Then.
O\(G) = A
r
JaG
Ot (~ )oii(~ ( ~ )) op(x(n )dl an
on
+ o(A).
Let G be a class of admissible sets, The domain G is said to be admissible . i. e,. G e G if the following conditions hold: (C-3 ) G is a bounded open set in R(2 ) such that D1 cGcD 2 : (C-4) G is the interior region of a C2-parameterized curve given by {x( ~ ) I ~ e [0. 1] } , With the definition of a class of admissible
where p(x ) is the solution to (vp ( x )) ' vcp(x )dx = -H[cp] , {~ G
~
(15)
i
- H[ii(G ) ] } .
( 16 )
for any cp e H6 (G). dl~ denotes the line element and ll..:.2 ~ on
v(' )
(18)
Spatial Doma in Identifi cation of Distribut ed Parame t e r Sys t ems
We shall implement an algorithm for finding - (G* )=0. As stated the domain G* such that 8~ in Lemma . 8~ (G) is obtained only through the boundary curve x ( e ) of the domain G. Then. for l>O. set as _ aJ (G ) xn+1 (e )=xn (e)-l 'Qa n n (e ) ( 19a ) ( 19b ) where Gn is the domain obtained by the n-th recursive computation and the boundary of Gn is denoted by aGn={xn (e ) leE [0 . 1] }. In (19b ). the operator PG is a projection on G. Denoting the boundary integral on aG by <• . • >oG . i . e . . for any ~ and ~.
we may write 8J .(G.OI ) =
aj .( G)
-~p ( x ) = -h(x ) ' { ~(t)-H[u ( Gt ) ] } for x E Gt with the boundary condition for
e E [ 0. 1]
( 25a) ( 25b )
where u(Gt ) is the solution to ( 24 ) and h( x ) =[h (x ) . h (x ) . . • • . hm (x )] ·. By solving 1 2 ( 24 ) and ( 25 ) . a pair of the solutions {u (Gt ). p(G t ) } is obtained . We (Step-3 ) Es timat ed domain co ns true tion shall set the perturbation parameter l in ( 23 ) as l=k . Then. the perturbation formula ( 23 ) is rewritten as
-
-
xt+k (e)=xt (e )-k
aut ~t (e ) ) an
( 26a) ( 26b)
< 01
•
aj .( G) --aa-
( 21)
>aG
Thus . comparing ( 21 ) with ( 15 ) . we obtain
--aa-
651
a ~ ( x ( e ))
on
op( x(e )) on
(22)
Hence (19 ) is _ o~ ( xn ( e )) op ( xn (e )) xn +1 (e )=x n (e)-l an an n(e ) ( 23a) ( 23b ) Based on ( 23 ). the following procedure to perform_computer implementation is stated : Let Gt be the estimated domain of Goat time t. We denote the boundary of Gt by aGt= { ~t ( e ) I eE[O.l] }. (Step-I ) Preassignment of an initially estimated domain At time to we select an initially estimated curve ~t ( e ) E G such that
o
~t ( O )· ~t ( 1 ) and construct the estimated doo _ 0 main Gt with the boundary oG t = { ~t ( e ) leE
o
0
0
[0 . 1]}. Let t be the discrete-time index t=tO+ik . i=O . 1. 2.3 . • . • • • . where k is the partitioned time . By using a newly observed data at each time t . recursive calculations for computing the estimated boundary curve associated with the domain Gt ( t = to . to+k. t o+2k . • • • ) are performed along Step-2 and Step-3.
The estimated domain Gt +k is determined by the curve oGt+k= { ~t+k ( e ) leE [0 . 1] } . Convergence properties of Gt are examined along the theoritical aspect which will be stated in the following section . Remark: The projector PG in ( 26b ) implies that the boundary curve ~t+k ( e ) is given by joining together the curve ;t+k (e ) and the boundary curve of D . Thus . it may reasonably be as2 sumed that the boundary curve ~t ( e ) is modified so that ;t (e ) may become C2-curve . 4.
CONVERGENCE PROPERTIES OF ESTIMATION ALGORITHM
It is difficult to examine asymptotic behaviors of the estimated domain Gt because the domain estimate i s performed by the estimated boundary curve { ~ ( e ) leE [0 . 1] } based on the solutions to ( 24 ) and ( 25 ). A possible method for investigating asymptotic behaviors is an evaluation of the cost functional jt (Gt ). Starting with ( 8 ) and ( 11 ) and noting that . from (14 ) . jt+k (Gt+k )-jt (Gt +k ) ( 29 ) 8jt (Gt )
jt+k (Gt+k )-jt (Gt ) +
we obtain
(Step - 2 ) Steady s tate solution of th e sy s t em and its adjoint one . Using the estimated domain Gt for G in ( 9 ). we write -.:.iI( x ) = f (x ) for x E Gt (t=tO+ik ) ( 24a ) with the boundary condition for
e E [ 0.1]
( 24b )
On the other hand . the adjoint equation ( 16 ) is expressed by
From Lemma and ( 20 ) . it follows that 8J (G ) = -k(- ( O~ ( x ( e )) ap(x(f ))) 2dl t t JaG on an t
+
o(k) .
e
( 31 )
Y. Sunahara, Sh . Aihara and F. Kojima
652
L W(t . q )
Then. from ( 30 ) and ( 31).mu1tiplying both sides of ( 30) by 1/ k. integrating the result over [to ' T] and letting k-+ O. we have -8 (t )dt + ~ H[u(G t )] '{ ~ y(t) - H[u(t . Go )] }dt - ~ H[u(Gt)] ' dv(t ) (32 )
dJ(t)
with
xexp{
r_
J,oG
( ou(~(e» on
op(x(f) ) )2dl
on
From (P-l) to (P-3 ). it follows that
(34 )
e
t
The following proposition states ~n asymptotic behavior of the solution pr~cess J (t ) of ( 32 ): [Proposition] The solution J (t ) to ( 32 ) satisfies the following equation. i,e , . for any ~>O . lim P[ 1 J (t ) t-+ 00
Proof
+
~ y(t) ' y(t ) 1 2> d 2t
Set as
q(t) = J (t ) + ~ y(t) ' y(t ) , 2t
= 0, (35) (36 )
<0
-28 (t)q (t )
if q(t).;, 0
-8(t )dt+y(t )dt+ry(t )' dv (t)
L
( 42 )
W(t .q ) < 0 .
( 44 )
See Has ' minskii
Consequently. for any ~ > O . 1980 . for more details ) lim P[ t-+ 00
1
q (t )
2 >~ ] = 0 .
1
( 45 )
Thus . the proof has been completed . ~onvergence property of the estimated domain Gt to the true domain Go is shown by the following theorem .
~ > O,
[Theorem-3] For any holds :
Using Lemma. the q (t )-process satisfies dq(t )
(41 )
-I y( t ) 1 (q (t) )2+ 2y (t )q(t )- 1 y( t ) 1 = - 1 y(t ) 1 {( q(t ) )2_ ,2~ i ~ l , q (tH } :o; O. ( 43 ) Hence . if q (t)';' 0
where O(t)=
f,oo;. (s) 1 ds }. t
and
!J
{-28(t)q(t) + 2y (t)q (t )
:0;
- 1 y ( t ) 1 (q ( t) ) 2_ 1 y ( t ) 1 }
lim P[ t-+ 00
( 37a )
the following limit
-
11
H[u(Gt )] - H[u (Go )]
11
(m»d=O.
R
(46 )
with -
1
'
J (t O) + - 2 y(t O) y(t O) . ( 37b ) 2t
Proof :Noting that lim v~t ) = 0 t-+ 00
where y (t ) and ry (t ) are a scalar and an m-dimensional process respectively given by y(t) = -~ [~y ( t ) - H[U(Gt ) ]] ' [~ y(t ) and
- H[u(t;G )]] + ~2 o 2t
( 39 )
lim t-+oo
~
lim t-+ 00
q(t)
0 ;
(P-2)
if q(t)=O . then 8 (t )=0. y(t )= ~ 2t2 and ry(t)=O
-
<
00
= 1
]
and
to
P{
L ry (t) , ry (t )dt < 00
t-+
00
=
Associated with (P-3). define W(t.q ) = { (q(t) )2 + LOC>I y(s )
t
oo
1
ds
t
ry (s) ' ry(s )ds }exp{Loo, t
00
respectively . From ( 35 ). ( 49 ) and ( 50 ). the equality ( 46 ) holds . DIGITAL SIMULATION STUDIES
Exampl e -I : The distributed system input f (x )
to
+L
t ........
+ ~ H[u(Gt )]'H[u(G t )] } w, p. l (50)
5. J
00
-
lim J (t ) = lim { -H[u(Gt )]'H[u(G o )]
if q(t)';' O. then 8(t » 0
P[ LOC>I y ( t) 1 dt
( 49)
w. p. 1
and
(P-l)
(P-3)
j(tH[u (s;G )]ds 0 0
it is easy to show that
With (C-1) and (C-2 ) . it can be shown that the q (t)-process has the following properties : ~
( 47 )
and
( 38 )
ry (t ) = ~ { ~ y(t) - H[u(Gt )] }.
w. p.1
y (s ) 1
ds } . ( 40)
and write the differential generator by L associated with (37 ) , we obtain
in (1) is given by f(x ) {=-2 forxEDicD1 ( 51a ) = 0 for x E D1 and -2 :0; f (x ) :0; 0 for xED 1 - Di ' (51b ) Letting m=l . i . e .• Wl=D l in ( 2). observations are made in the form of dy (t ) =
(j[ h1 (x )u(t.x )dx)dt + dv (t ) D1
(52)
653
Spa tial DomaiR Identificati on of Distributed Pa rame t e r Systems
with y(O)=O . where for x e Di . h (x)f=h (53a) 1 Lo for x e Dl and Os; h 1 ( x ) s; h forxeD1-Di . ( 53b) (Preassignment of Go and D1 ) In Example-I . for the purpose of simulation experiments. the true spatial domain Go which is to be identified and the observation domain Dl were preassigned as shown in Fig . 1. The domain Go was decomposed by 101 triangles of subregions with 64 nodes . curve of true domain Go
X
la'
dye tl/ d t
100
50
o
50
100
Time
Fig . 2. Sample run of the scalar observation dy (t )/dt in Example-l 150
(52). respectively. obtain the solutions u and p to ( 24) and (25 ). Step-3 : By substituting the solutions and p obtained in Step-2 into ( 26 ). the k-th iterated curve of the estimated boundary ~t +k (E ) is computed . 0 St ep-4 Construct the estimated domain Gt +k
u
100
50
o
o
50
10 0
150
200
Xl
Fig . 1. Unknown domain Go ' given region Dl and their triangulations in Example-l (Observation data ) In simulation experiments. the observation data dy (t ) were generated by ( 52 ) and ( 53 ) . where the u( t . x )-process was obtained by applying the finite element method to ( 1 ) with the input given by ( 51 ) and the s tandard Brownian process v(t ) generated by a usual way. A sample run of the observation process is shown in Fig. 2. (Domain identificaion ) Based on the computer algorithm stated in the previous section . the numerical procedure is as follows : St ep-l Preas sign the initial curve { ~t (E)} as shown in Fig . 3. 0 St ep-2 With the input f and the observation data . i.e .. the output data. given by ( 51 ) and
based on the boundary curve derived in Step-3 . Letting k=I.2 . • • • . Steps ( 2 ) to ( 4) give a forwardly recurre~t algorithm to obtain the running estimate Gt with the input f ( x ) and the output dy(t ) as a set of the given data . Figure 3 shows the evolution of estimated boundary curves .
Example-2
The distributed input f (x l is
given by
f(x ) [ :
-40 0
for x e Di for x e Dl
and -40 S; f ( x) S; 0 for x e D1- Di ( 54b ) The observation mechanism is set with fiftyone observation domains. loll' 101 2 , •• 101 51 , (Preassignment of Go and Dl in Example-2 ) Regarding the true domain Go and the preassigned domain D1 . the same decompositions as in Example-I were adopted. These are illustrated in
1500
15 .0 iteration 1580 curve
,, ,,
lOO
,
··· ./
5.0
,,'/
,// ./
o Fig . 3.
( 54a )
15.0 5.0 10.0 20.0 Estimated boundary of unknown domain Go in Example-l
Y. Sunahara, Sh. Aihara and F. Kojima
654
Chavent . G ( 1978) . About the identification and modeling of miscible or immisicible displacements in porous media In Lec ture Not es in Control and Information Sci ences, Vol 1 . pp . 196-220 . Springer. New York . Chenais . D (1975 ) On the existence of a solution in a domain identification problem j. Math. Anal. Appl. 52. pp 189-219 Delfour . M.. G Payre and J-P Zolesio ( 1983 ) Optimal design of a minimum weight thermal diffuser with constraint on the output ther mal power flux . Appl. Math. Optim. 9 . pp , 225-262 St ochas ti c StabiliHas' minskii . R. Z. ( 1980) . ty of Differ ential Equations , Sijthoff and Noordhoff . Alphen aan den Rijn Pironneau . 0 (1977 ) Variational methods for the numerical solutions of free boundary problems and optimum design problems . In Control Th eory of Systems Governed by Partial Differ ential Equations, 209 . Academic Press . New York . Pironneau . 0 (1982 ). Optimal shape design for elliptic systems . In Lec ture Not es in Control and Information Sci ences , Vo1. 38 . pp , 42-66 . Springer New York . ACKNOWLEDGMENTS The authors wish to e xpress their thanks to Messrs . Teruo OKA and Toshiyuki MORIOKA for their helpful assistance .
Fig
4 . Each observation domain \'\ (i=1.2.· . 51 ) was given by a triangle as shown in Fig 5 Using the 51-dimensional observation mechanism . a sample run of the estimated domain is shown in Fig 5. iterated by the same procedure as described in Example-I . CONCLUSIONS A feasible method has been presented in this paper. for identifying the boundary curve of the unknown domain of distributed parameter systems under noisy observations The following three aspects in this paper should be emphasized as concluding remarks ( 1 ) Math ematical model Recognizing the fact that dynamic behaviors of many real physical systems are . in fact. distributed. a mathematical model is described by partial differential equations of parabolic type not elliptic . Furthermore. noisy observations are taken into account from viewpoints of practical applications . (2 ) Recursive algorithm for computing th e unknown spatial domain By using current information of the input-output data . the unknown spatial domain can be estimated by recursive computations . (3 ) Us e of th e stationary state information of th e system A difficulty arises in computing recursively the unknown spatial domain because of increase of the observation data . Consequently . use is made of the steady state solution of the system for the purpose of circumventing massive computations . Studies on convergency properties described in Section 4 support the feasibity and the validity of the proposed method . Turning our eyes to context of shape design by Pironneau (1977 . 1982 ) and Begis . et al ( 1975 ). it can be understood that the idea involved in this paper bridges between the method of the optimal shape design and that of maximum likelihood estimate
domain Go region Dl
50
REFERENCES Begis . D. and R Glowinski (1975 ). Application de la method des e lements fini a l ' approximation d ' un probleme de domain optimal . Methodes de resolution des problemes approches Appl. Math . Opti m. 2 . pp . 130 - 169 .
o
50
100
150
iteration 7500 15 .0
true curve
,
:' 10.0 :,
,,
5 .0
'-'
o
5.0
Fig. 5.
100
200
X,
Fig . 4. Unknown domain Go ' given region D1 and their triangulations in Example-2
15.0
200
Xl
Estimated boundary of unknown domain Go in Example-2