A relaxation scheme for Hamilton–Jacobi–Bellman equations

Applied Mathematics and Computation 186 (2007) 806–813 www.elsevier.com/locate/amc

A relaxation scheme for Hamilton–Jacobi–Bellman equations

q

Shuzi Zhou *, Zhanyong Zou Institute of Applied Mathematics, Hunan University, Changsha, China

Abstract In this paper we propose a relaxation scheme for solving discrete Hamilton–Jacobi–Bellman equations based on Scheme I of Lions and Mercier. The convergence of the new scheme has been established. Numerical example shows the scheme is efficient. The existence of the solution of Hamilton–Jacobi–Bellman equation has also been discussed.  2006 Elsevier Inc. All rights reserved. Keywords: Relaxation scheme; Hamilton–Jacobi–Bellman equation; Convergence; Existence

1. Introduction We consider the following discrete Hamilton–Jacobi–Bellman equation: find u 2 Rn such that max fAj u  f j g ¼ 0;

16j6m

ð1:1Þ

where Aj 2 Rn·n, f j 2 Rn, j = 1, . . . , m. Eq. (1.1) is obtained by discretizing continuous Hamilton–Jacobi– Bellman equation. See [1,2] and the references therein. Lions and Mercier [1] has proposed two numerical schemes for solving (1.1). See also [3]. Scheme I Step 1. Given e > 0, k:¼1, for some j find u0,m such that Aj u0;m ¼ f j : Step 2. Let N = (k  1) m,uN,0 = uN,m. For j = 1, . . . , m find uN+j,j such that maxfAj uN þj;j  f j ; uN þj;j  uN þj1;j1 g ¼ 0: Step 3. If kukm,m  uN,0k < e then output ukm,m otherwise k :¼ k + 1, go to Step 2. q *

Supported by NNSF of China (No. 10571046). Corresponding author. E-mail address: [email protected] (S. Zhou).

0096-3003/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.08.025

S. Zhou, Z. Zou / Applied Mathematics and Computation 186 (2007) 806–813

807

Assume Aj ¼ ðajls Þ; f j ¼ ðflj Þ. Let p

Aðp1 ; . . . ; pn Þ ¼ ðalsl Þ;

p

f ðp1 ; . . . ; pn Þ ¼ ðfl l Þ:

ð1:2Þ pl

That is: the lth row of matrix A(p1, . . . , pn) is the lth row of matrix A , the lth component of vector f(p1, . . . , pn) is the lth component of vector f pl . Now we formulate Scheme II of Lions and Mercier in the notation above. Scheme II Step 1. k :¼ 0, for some j find u0 such that Aj u 0 ¼ f j :

ð1:3Þ pkl

Step 2. For l = 1, . . . , n find such that   j k j k k j j pl ¼ min j 2 f1; . . . ; mg : ðA u  f Þl ¼ max fðA u  f Þl g : 16j6m

ð1:4Þ

Step 3. Compute uk+1 as the solution of Aðpk1 ; . . . ; pkn Þukþ1 ¼ f ðpk1 ; . . . ; pkn Þ:

ð1:5Þ

Step 4. If uk+1 = uk then output uk otherwise k:¼k + 1 and go to Step 2. In the last decades many numerical schemes have been given for solving (1.1). But the above schemes are still playing very important role. See [3–6] and the references therein. In this paper we propose, based on Scheme I above, a relaxation scheme with a parameter x, which for x = 1 is just Scheme I. In our numerical example, the new scheme with x = 0.5, 0.8, 0.9 is faster than Scheme I (x = 1). The monotone convergence of the new scheme has been proved. On the other hand, we give an existence theorem for (1.1), the proof of which is based on Scheme II. 2. Assumptions and existence of solution In [6] we proposed the following conditions for (1.1). Condition A: Aj, j = 1, . . . , m, are L-matrixes and have strict diagonal dominance (see [7]). Condition B: Aj, j = 1, . . ., m, are L-matrixes, and all the matrixes A(p1, . . . , pn), pl = 1, . . . , m, l = 1, . . . , n, are irreducible, weakly diagonally dominant. In [6] we have proved the following theorem. Theorem 2.1. If Condition A or Condition B holds then (1.1) has a unique solution. In this paper we consider weaker condition as follows. Condition A*: All the matrixes A(p1, . . . , pn), pl = 1, . . . , m, l = 1, . . . , n, are M-matrixes. The following lemma is obvious. Lemma 2.1. If Condition A or Condition B holds then Condition A* holds. Example 2.1. n = k = 2, and    3 1 5=2 A1 ¼ ; A2 ¼ 4 2 1

1 1

 :

It is easy to see Condition A* holds. But A1 is not diagonally dominant. Hence neither Condition A nor Condition B holds. Now we prove that the weaker Condition A* still implies the existence of solution of (1.1). Theorem 2.2. If Condition A* holds then (1.1) has a unique solution.

808

S. Zhou, Z. Zou / Applied Mathematics and Computation 186 (2007) 806–813

Proof 1. Consider Scheme II under Condition A*. Since all A(p1, . . . , pn), pl = 1, . . . , m, l = 1, . . . , n, are Mmatrixes, uk, k = 0, 1, . . ., in Scheme II are well defined, First, we prove uk is decreasing monotonically, i.e.,    6 ukþ1 6 uk 6    6 u1 6 u0 :

ð2:1Þ

By (1.5) we have Aðp01 ; . . . ; p0n Þu1 ¼ f ðp01 ; . . . ; p0n Þ;

ð2:2Þ

which combining with (1.3) and (1.4) yields Aðp01 ; . . . ; p0n Þu0  f ðp01 ; . . . ; p0n Þ P Aj u0  f j ¼ 0 ¼ Aðp01 ; . . . ; p0n Þu1  f ðp01 ; . . . ; p0n Þ: 1

ð2:3Þ

0

Since Aðp01 ; . . . ; p0n Þ is an M-matrix, (2.3) means u 6 u . Similarly, we have by (1.5) that Aðp11 ; . . . ; p1n Þu2 ¼ f ðp11 ; . . . ; p1n Þ; which combining with (1.4) and (2.2) yields Aðp11 ; . . . ; p1n Þu1  f ðp11 ; . . . ; p1n Þ P Aðp01 ; . . . ; p0n Þu1  f ðp01 ; . . . ; p0n Þ ¼ Aðp11 ; . . . ; p1n Þu2  f ðp11 ; . . . ; p1n Þ: Hence we have u2 6 u1. Then we derive (2.1) by induction. It follows from (1.4) and (1.5) that max fAj uk  f j g ¼ Aðpk1 ; . . . ; pkn Þuk  f ðpk1 ; . . . ; pkn Þ ¼ Aðpk1 ; . . . ; pkn Þðuk  ukþ1 Þ;

16j6m

k ¼ 0; 1; . . .

ð2:4Þ

Since the set {(p1, . . . , pn): pl = 1, . . . , m, l = 1, . . . , n} is a finite set there exist positive integers q and k with q > k such that ðpq1 ; . . . ; pqn Þ ¼ ðpk1 ; . . . ; pkn Þ: Therefore, we have Aðpq1 ; . . . ; pqn Þ ¼ Aðpk1 ; . . . ; pkn Þ; f ðpq1 ; . . . ; pqn Þ ¼ f ðpk1 ; . . . ; pkn Þ: Then by (1.5) we obtain uqþ1 ¼ ukþ1 ; which and (2.1) results in uqþ1 ¼ uq ¼    ¼ ukþ2 ¼ ukþ1 :

ð2:5Þ

It follows from (2.4) and (2.5) that max fAj ukþ1  f j g ¼ 0;

16j6m

which means uk+1 is a solution of (1.1). The existence of solution has been proved. Finally, we prove the uniqueness of solution. Assume u and u* are solutions of (1.1), i.e., max fAj u  f j g ¼ 0;

ð2:6Þ

max fAj u  f j g ¼ 0:

ð2:7Þ

16j6m 16j6m

It is easy to see from (2.6) and (2.7) that there exist (p1, . . . , pn) and ðp1 ; . . . ; pn Þ such that Aðp1 ; . . . ; pn Þu  f ðp1 ; . . . ; pn Þ ¼ 0;

ð2:8Þ

Aðp1 ; . . . ; pn Þu  f ðp1 ; . . . ; pn Þ ¼ 0; Aðp1 ; . . . ; pn Þu  f ðp1 ; . . . ; pn Þ 6 0; Aðp1 ; . . . ; pn Þu  f ðp1 ; . . . ; pn Þ 6 0:

ð2:9Þ ð2:10Þ ð2:11Þ

(2.8) and (2.11) implies u* 6 u. But (2.9) and (2.10) implies u* P u. Hence u* = u. The proof is complete.

h

S. Zhou, Z. Zou / Applied Mathematics and Computation 186 (2007) 806–813

809

Remark 2.1. (2.4) means that if uk+1 = uk then uk is the solution of (1.1). This shows the reasonability of Step 4 of Scheme II. Now we consider another condition weaker than Condition B. Condition C: Aj, j = 1, . . . , m, are irreducible, diagonally dominant M-matrixes. The following example shows that Condition C cannot guarantee the existence of solution. Example 2.2. n = k = 2, and    2 1 1 1 2 A ¼ ; A ¼ 1 1 1 1

1 2

 ;

T

2

f ¼ f ¼ ð1; 2Þ : Obviously, Condition C holds in this example. On the other hand, Condition A* does not hold since A(2, 1) is not an M-matrix but a singular one. We prove in this case (1.1) has no solution. If it is not true then there T exists u ¼ ðu1 ; u2 Þ 2 R2 satisfying maxfAj u  f j g ¼ 0;

ð2:12Þ

j¼1;2

i.e. maxf2u1  u2  1; u1  u2  1g ¼ 0;

ð2:13Þ

maxfu1

ð2:14Þ

þ

u2

þ

2; u1

þ

2u2

þ 2g ¼ 0:

Then u* must satisfy one of the following four equations: Aðp1 ; p2 Þu  f ðp1 ; p2 Þ ¼ 0;

p1 ; p2 ¼ 1; 2;

i.e. 

2 1  2 1  1 1  1 1

   1 1  ¼ 0; u  2 1    1 1  ¼ 0; u  2 2    1 1  ¼ 0; u  2 1    1  1 u  ¼ 0; 2 2

ð2:15Þ ð2:16Þ ð2:17Þ ð2:18Þ

The solution of (2.15) is (1, 3)T. But it does not satisfy (2.13). Both (2.16) and (2.18) has the same solution (0, 1)T. But it does not satisfy (2.14). (2.17) has no solution at all. Therefore, u* is not a solution of (2.12). This contradiction shows that (2.12) has no solution. 3. Relaxation scheme and convergence We propose a relaxation scheme which is an extension of Scheme I. Scheme RS Step 1. Given e > 0, x 2 (0, 1], k :¼ 1, for some j find u0,m such that Aj u0;m ¼ f j :

ð3:1Þ N,0

N, m

Step 2. Let N = (k  1)m, u = u , j :¼ 1. Step 3. Compute vN+j,j as the solution of maxfAj vN þj;j  f j ; vN þj;j  uN þj1;j1 g ¼ 0:

ð3:2Þ

810

S. Zhou, Z. Zou / Applied Mathematics and Computation 186 (2007) 806–813

Step 4. Compute uN þj;j ¼ ð1  xÞuN þj1;j1 þ xvN þj;j :

ð3:3Þ

If j = m then go to Step 5 otherwise j :¼ j + 1 and go to Step 3. Step 5. If kukm,m  uN,0k < e then output ukm,m otherwise k :¼ k + 1 and go to Step 2. We have the following convergence theorem. Theorem 3.1. Assume that Condition A* holds, and that ukm,m, k = 0, 1, 2, . . ., are produced by Scheme RS. Then ukm,m is monotonely decreasing and convergent to the solution of (1.1). Proof 2. It follows from (3.2) that for j = 1, . . . , m holds vðk1Þmþj;j 6 uðk1Þmþj1;j1 ;

k ¼ 1; 2; . . .

ð3:4Þ

Letting k = 1 in (3.4) we obtain vj;j 6 uj1;j1 ;

j ¼ 1; . . . ; m;

ð3:5Þ

which and (3.3) implies that uj;j 6 uj1;j1 ;

j ¼ 1; . . . ; m;

i.e. um;m 6 um1;m1 6    6 u0;0 ¼ u0;m :

ð3:6Þ

By similar argument we derive u2m;m 6 u2m1;m1 6    6 um;0 ¼ um;m :

ð3:7Þ

It follows from (3.6), (3.7) and induction that for k = 1, 2, . . . holds ukm;m 6 ukm1;m1 6    6 uðk1Þm;0 ¼ uðk1Þm;m 6 uðk1Þm1;m1 6    6 uðk2Þm;0 :

ð3:8Þ

It is well known that linear complementarity problem (3.2) is equivalent to a obstacle problem and uN+j1,j1 is the !0 obstacle !1 #, and that the solution of obstacle problem is monotone respect to the obstacle (ref. Lemma 4.2 in [6]). Hence (3.2) and (3.8) implies vkm;m 6 vðk1Þm;m ; vkm1;m1 6 vðk1Þm1;m1 ; . . . ; vðk1Þmþ1;1 6 vðk2Þmþ1;1 : Now we prove {u j

km,m

} has a bound below. Since Condition

A*

ð3:9Þ

holds, (1.1) has a unique solution ~u, i.e.

j

max fA ~ u  f g ¼ 0:

16j6m

ð3:10Þ

It follows from (3.1) that max fAj u0;m  f j g P 0:

16j6m

Hence there exist indexes p1, . . . , pn such that Aðp1 ; . . . ; pn Þu0;m  f ðp1 ; . . . ; pn Þ P 0:

ð3:11Þ

On the other hand, (3.10) implies that u  f ðp1 ; . . . ; pn Þ 6 0: Aðp1 ; . . . ; pn Þ~

ð3:12Þ

By Condition A*, we know A(p1, . . . , pn) is an M matrix. Therefore, we have by (3.11) and (3.12) that u: u0;m P ~

ð3:13Þ

Now we prove u: v1;1 P ~

ð3:14Þ

S. Zhou, Z. Zou / Applied Mathematics and Computation 186 (2007) 806–813

811

e ¼ f1; . . . ; ng; I ¼ fs 2 N e : v1;1 ¼ u0;m g. It follows from (3.13) that Let N s s 0;m P~ us v1;1 s ¼ us

for

s 2 I:

ð3:15Þ

By (3.2) we have ðA1 v1;1  f 1 Þs ¼ 0

for

e n I: s2N

ð3:16Þ

On the other hand, (3.10) implies A1 ~ u  f 1 6 0:

ð3:17Þ

Denote by fP the subvector of f corresponding to the index set P, and by AP,Q the submatrix of A corresponding to the row and column index set P and Q. Then (3.16) can be rewritten as 1;1 1 1 A1~ e v1;1 ~ nI;I vI ¼ fN ~ nI : ~ þ AN N nI; N nI N nI

ð3:18Þ

By (3.17) we have A1~ e ~ uN~ nI þ A1N~ nI;I ~ uI 6 fN~1nI ; N nI; N nI which combining with (3.18) yields A1~ e ðv1;1  ~ uÞN~ nI þ A1N~ nI;I ðv1;1  ~ uÞI P 0: N nI; N nI

ð3:19Þ

Since A1N~ nI;I 6 0; AN~ nI;N~ nI is an M-matrix, we have by (3.15) that ðv1;1  ~ uÞN~ nI P 0:

ð3:20Þ

(3.15) together with (3.20) is just (3.14). By (3.14), (3.3) and x 2 (0, 1] we derive that u: u1;1 P ~ It is easy to prove by similar way and induction that vkmþj;j P ~ u; ukmþj;j P ~ u;

j ¼ 1; . . . ; m;

k ¼ 0; 1; 2; . . . *,j

ð3:21Þ *,j

It follows from (3.8), (3.9) and (3.21) that these exist v , u , j = 1, . . . , m, such that vkmþj;j ! v;j ;

ukmþj;j ! u;j ; ðk ! 1Þ;

ð3:22Þ

and v;j P ~ u;

u;j P ~ u:

ð3:23Þ

Letting k ! 1 in (3.8) we obtain u;m 6 u;m1 6    6 u;0 ¼ u;m : Hence we have u;1 ¼ u;2 ¼    ¼ u;m ; which and (3.3) implies v;1 ¼ v;2 ¼    ¼ v;m ¼ u;1 ¼    ¼ u;m :

ð3:24Þ

By (3.2) and (3.24) we derive Aj v;m  f j 6 0;

j ¼ 1; . . . ; m:

Hence we have max fAj v;m  f j g 6 0:

16j6m

ð3:25Þ

812

S. Zhou, Z. Zou / Applied Mathematics and Computation 186 (2007) 806–813

Using the similar way for proving (3.14), we derive from (3.25) and (3.10) that u; v;m 6 ~ which combining with (3.23) and (3.24) yields u; u;m ¼ ~ i.e. u; ukm;m ! ~

ðk ! 1Þ:

The theorem has been proved.

h

Corollary 1. Assume Condition A or Condition C holds. Then we have the conclusion same as in Theorem 3.1. 4. Numerical example We use example 1 in [3], i.e. n = m = 2, X = (0, 1) · (0, 1), A1 ¼ 

o2 o2 o2   0:5 ; ox2 ox oy oy 2

o2 o2 o2   0:1 ; ox2 ox oy oy 2 u ¼ xð1  xÞyð1  yÞ; A2 ¼ 0:5

f 1 ¼ f 2 ¼ maxfA1 u; A2 ug: The discretization of the above second order derivatives are: o2   h2 Dþ h;x Dh;x ; ox2

o2   h2 D h;y d h;y ; oy 2

o2 1 þ    h2 ½Dþ h;x Dh;y þ Dh;x Dh;y ; ox oy 2  where D h;x ; Dh;y denote the forward and backward difference respectively in x and y, h = 1/10, 1/20. We use Algorithm RS to solve the discrete problem. The subproblems (3.2) are solved by Scheme I. Take e = 105, x = 0.1, 0.5, 0.8, 0.9, 1.0 and 1.1, respectively. Table 1 shows the 1-norm of the residual R = max16j6m{Ajuk,m  f j} when iteration terminates. We see that R  0 for x 6 1 and R is big for x = 1.1. It appears also for x = 1.3, 1.5, 1.8, 1.9. Table 2 shows the relation between iteration number k and relaxation number x. Tables 3 and 4 show the value of uk,m at (x, y)T = (0.5, 0.5)T for h = 1/10 and h = 1/20 respectively.

Table 1 x

0.1

0.5

0.8

0.9

1.0

1.1

kRk1 1 h ¼ 10 1 h ¼ 20

1.637 · 105 3.244 · 104

2.426 · 1014 3.390 · 1013

3.606 · 1014 3.988 · 1013

2.897 · 1014 2.897 · 1014

2.349 · 1014 3.382 · 1013

1.853 · 102 2.786 · 102

Table 2 x

0.1

0.5

0.8

0.9

1.0

k 1 h ¼ 10 1 h ¼ 20

200 300

85 200

50 170

47 201

169 281

S. Zhou, Z. Zou / Applied Mathematics and Computation 186 (2007) 806–813

813

Table 3 0.1

x k,m

0.5

0.8

0.9

1.0

0.074597215 0.069596631 0.064864841 0.064381072

0.069237391 0.065471011 0.064686484 0.064381072

0.067450783 0.064944043 0.064527151 0.064381072

0.065664174 0.065655725 0.065601722 0.064381072

0.5

0.8

0.9

1.0

0.072935081 0.066534959 0.064284738 0.064214317

0.067312343 0.064476811 0.064283484 0.064214317

0.067475812 0.064605112 0.064282158 0.064214317

0.065449468 0.065215940 0.064284998 0.064214318

T

u at (0.5,0.5) , h = 1/10 k=1 0.081743648 k=2 0.080074193 k=3 0.078530160 Last k 0.064381093

Table 4 x uk,m at (0.5,0.5)T, k=1 k=2 k=3 Last k

0.1 h = 1/20 0.081123414 0.079272722 0.077843915 0.064214323

We can see from Table 2 that the algorithm for x = 0.5,0.8,0.9 is faster that for x = 1.0. Table 3 displays the monotonicity of the algorithm. References [1] P.L. Lions, B. Mercier, Approximation numerique des equations de Hamilton–Jacobi–Bellman, R.A.I.R.O., Anal. Numer. 14 (1980) 369–393. [2] M. Boulbrachene, M. Haiour, The finite element approximation of Hamilton–Jacobi–Bellman equations, Comput. Math. Appl. 41 (2001) 993–1007. [3] R.H.W. Hoppe, Multigrid method for Hamilton–Jacobi–Bellman equations, Numer. Math. 49 (1986) 239–254. [4] M. Sun, Domain decomposition method for solving Hamilton–Jacobi–Bellman equations, Numer. Funct. Anal. Optim. 14 (1993) 145– 166. [5] M. Sun, Alternating direction algorithms for solving Hamilton–Jacobi–Bellman equations, Appl. Math. Optim. 34 (1996) 267–277. [6] S.Z. Zhou, W.P. Zhan, A new domain decomposition method for an HJB equation, J. Comput. Appl. Math. 159 (2003) 195–204. [7] D. Young, Iterative Solution of Large Linear Systems, Academic Press, New York, 1971.