View of On the convergence of a method for solving two point boundary value problems by optimal control

Rev. Anal. Num´er. Th´eor. Approx., vol. 31 (2002) no. 2, pp. 217–227 ictp.acad.ro/jnaat

ON THE CONVERGENCE OF A METHOD

FOR SOLVING TWO POINT BOUNDARY VALUE PROBLEMS BY OPTIMAL CONTROL

ERNEST SCHEIBER^∗

Abstract. Using the idea of the least squares method, a nonlinear two point boundary value problem is transformed into an optimal control problem. For solving the optimal control problem it is used the gradient method. The conver- gence of the method is investigated and numerical results are reported.

MSC 2000. 65P40.

Keywords. Two point boundary value problem, optimal control, least squares method, gradient method.

1. INTRODUCTION

In this paper we study the convergence property of a method to solve the nonlinear two point boundary value problem (NTPBVP)

x^(m)(t) =f x(t),x(t), . . . , x˙ ^(m−1)(t), t, t∈[a, b], (1)

m

X

j=1

αi,jx^(j−1)(a) +βi,jx^(j−1)(b)=γi, i∈ {1,2, . . . , m}

(2)

using an optimal control problem (OCP). For the problem x(t)¨ = f(x(t),x(t), t),˙ t∈[0, T], (3)

x(0) = α, (4)

x(T) = β, (5)

the method was described in a previous note [11].

Sokolowski, Matsumura and Sakawa [12] used optimal control methods to solve two point boundary value problems of the form

−_dt^dha t, y(t),^dy_dt^dy_dtⁱ+qy(t) =f(t), t∈[0,1], y(0) =y(1) = 0.

The nonlinear two point boundary value problems and the optimal control problems are connected. The necessary optimality conditions, as Pontryagin’s

∗“Transilvania” University of Bra¸sov, Faculty of Mathematics and Computer Science, Department of Computer Science, st. Iuliu Maniu 50, 2200 Bra¸sov, Romania, e-mail:

[email protected].

maximum principle, lead for some optimal control problem to a nonlinear two point boundary value problem such as (3)–(5).

The multiple shooting method (Keller H. B. [6], Marzulli P., [8]), the col- location method (Ascher U., Christiansen J., Russell R. D., [1], [2]) are well known and widely used to solve a NTPBVP.

In our case, the derived OCP may be solved efficiently using the gradient method. The application of the gradient method to solve optimal control problems is well known: Polak E. [10], Polak E., Klessig R., 1973; Fedorenko P. R., 1878 and Miele A. [9].

Another possible method to solve the optimal control problem is the control parametrization (Goh C. Z., Teo K. L., 1988, Teo K. L., Goh C. J., Wong K. H., 1991).

Although the NTPBVP (1)-(2) has not a very general form, thanks to the boundary conditions, our approach emphasizes a class of NTPBVP which may be efficiently solved using optimization techniques.

2. STATEMENT OF THE PROBLEM

Consider the NTPBVP (1)-(2).

We assume that the NTPBVP has an unique solution and that f is con- tinuous together with his partial derivates of first and second order. Ifx(t) is the solution of the NTPBVP (1)-(2) then the pair (u(t), x(t)) is the solution of the following OCP

(6) minimize I(u) = Z b

a

u(t)−f(x(t),x(t), . . . , x˙ ^(m−1)(t), t)²dt subject to

(7) x^(m)(t) =u(t), t∈[a, b],

and (2).

Denoting x1(t) = x(t), x2(t) = ˙x(t), . . . xm(t) = ˙xm−1(t), u(t) = ˙xm(t) and x= (x₁, x₂, . . . , xm),the above problem may be written as an OCP for a first order differential system:

(8) minimize I(u) = Z b

a

u(t)−f(x₁(t), x₂(t), . . . , x_m(t), t)²dt subject to

x(t) =˙ Qx(t) +ξ_mu(t), t∈[a, b], (9)

Ax(a) +Bx(b) =γ, (10)

where

Q=

















0 1 0 . . . 0 0 0 1 . . . 0 ... ... ... ... 0 0 0 . . . 1 0 0 0 . . . 0















 ,

ξ_m = (0,0, . . . ,0,1)^t, A= (α_i,j)_1≤i,j≤m, B = (β_i,j)1≤i,j≤m,

γ = (γi)1≤i≤m.

For givenu the solution of the linear system (9) is

(11) x^u(t) =H(t)c+

Z t a

ϕ(t, s)u(s)ds, where

H(t) =















1 ^t−a_1! ^(t−a)_2! ² . . . ^(t−a)_(m−1)!^m−1 0 1 ^t−a_1! . . . ^(t−a)_(m−2)!^m−2

... ... ... ...

0 0 0 . . . 1













 ,

ϕ(t, s) =^(t−s)_(m−1)!^m−1, ^(t−s)_(m−2)!^m−2, . . . , 1^t, c= (c₀, c₁, . . . , cm−1)^t.

Using the shooting method, in order to satisfy the boundary condition (10), the vectorc is the solution of the following algebraical system

[A+BH(b)]c=γ−B Z b

a

ϕ(b, s)u(s)ds.

We suppose that the matrix R=A+BH(b) is not singular. It results that (12) x^u(t) =H(t)R⁻¹γ+

Z _b

a

K(t, s)u(s)ds, whereK(t, s) =ϕ₊(t, s)−H(t)R⁻¹Bϕ(b, s) and

ϕ+(t, s) =

(ϕ(t, s), if a≤s≤t≤b 0, if a≤t < s≤b.

To solve the OCP (8)–(10) by the gradient method, it requires to construct the sequence

uk+1=uk−µkI⁰(uk)

starting with a function u₀ ∈ L₂[a, b]. The descent parameter µ_k is usually computed as the solution of the one dimensional optimization problem

I(u_k−µ_kI⁰(u_k)) = min{I(u_k−µI⁰(u_k)) :µ≥0}.

If we denote L(x, u, t) = [u(t)−f(x1(t), x2(t), . . . , xm(t), t)]² then the Gˆa- teaux derivate of the cost functional is

I⁰(u)(δu) = Z b

a

h

Lx(x^u(t), u(t), t), δx(t)+Lu(x^u(t), u(t), t)δu(t)ⁱdt, where the functionsδxand δu satisfy the linear boundary value problem

δx(t) =˙ Q δx(t) +ξ_mδu(t), t∈[a, b], A δx(a) +B δx(b) = 0.

From (12) it results that

δx(t) = Z b

a

K(t, s)δu(s)ds and then

I⁰(u)(δu) =

= Z b

a

Z _b

a

DL_x x^u(t), u(t), t, K(t, s)^Edt+L_u x^u(s), u(s), s

δu(s)ds.

Hence the expression of the gradient becomes I⁰(u)(s) =

Z b a

D

Lx x^u(t), u(t), t, K(t, s)^Edt+Lu x^u(s), u(s), s. For the problem (3)–(5) the gradient of the cost functional may be computed by

(13) I⁰(u) =Lu(x^u₁, x^u₂, u, t)−p^u₂,

where p^u₁ and p^u₂ are the solutions of the following two point boundary value problem (the co-state system)

p˙₁ = L_x1(x^u₁, x^u₂, u, t), (14)

p˙2 = −p₁+Lx2(x^u₁, x^u₂, u, t), (15)

p2(0) = 0, (16)

p₂(T) = 0.

(17)

In this case, for the control function u the corresponding trajectory is given by

x^u₁(t) =α+^β−α_T t+ Z t

0

(t−s)u(s)ds−_T^t Z T

0

(T−s)u(s)ds, (18)

x^u₂(t) = ^β−α_T + Z t

0

u(s)ds−_T¹ Z T

0

(T−s)u(s)ds.

(19)

From (14)–(17) it follows that

p^u₂(t) = (20)

= Z t

0

hL_x2(x^u₁(s), x^u₂(s), u(s), s)−(t−s)L_x1(x^u₁(s), x^u₂(s), u(s), s)ⁱds−

−_T^t Z T

0

hL_x2(x^u₁(s), x^u₂(s), u(s), s)−(T−s)L_x1(x^u₁(s), x^u₂(s), u(s), s)ⁱds.

3. THE CONVERGENCE RESULT

We state a convergence result for the method considered above applied to the problem (3)–(5).

Ifu0∈L2[0, T] we denote by M_I(u0) the set defined by M_I(u0)={u∈L₂[0, T] :I(u)≤I(u₀)}

and we introduce the assumption:

(H) For any u, v∈M_I(u0) there existsC >0 such that

f(x^u₁, x^u₂, t)−f(x^v₁, x^v₂, t)≤Cku−vk₂;

u_∂x^∂f

k(x^u₁, x^u₂, t)−v_∂x^∂f

k(x^v₁, x^v₂, t)≤

≤C|u(t)−v(t)|+|x^u₁(t)−x^v₁(t)|+|x^u₂(t)−x^v₂(t)|;

f(x^u₁, x^u₂, t)_∂x^∂f

k(x^u₁, x^u₂, t)−f(x^v₁, x^v₂, t)_∂x^∂f

k(x^v₁, x^v₂, t)≤

≤C|x^u₁(t)−x^v₁(t)|+|x^u₂()−x^v₂(t)|

for anyt∈[0, T] and k∈ {1,2}.

If the functionf has continuous and bounded first and second order partial derivatives then the assumption (H) is satisfied.

We state some preliminary results.

Theorem 1. There exist the positive constants C₁ andC₂ such that

|x^u₁(t)−x^v₁(t)| ≤C₁ku−vk₂,

|x^u₂(t)−x^v₂(t)| ≤C₂ku−vk₂, for any u, v∈L2[0, T]and any t∈[0, T].

Proof. From (18) we find

x^u₁(t)−x^v₁(t) = Z t

0

(t−s)[u(s)−v(s)]ds− _T^t Z T

0

(T−s)[u(s)−v(s)]ds.

It follows that

|x^u₁(t)−x^v₁(t)| ≤ Z t

0

(t−s)²ds

¹₂ Z t 0

[u(s)−v(s)]²ds ¹₂

+ +

Z _T

0

(T−s)²ds

¹₂ Z _T

0

[u(s)−v(s)]²ds ¹₂

≤²

√ 3

3 T³²ku−vk₂. Thus C1 = ²

√3

3 T³²ku−vk₂. Analogously, we deduce

|x^u₂(t)−x^v₂(t)| ≤C2ku−vk₂, withC2 = 1 +

√3 3

T¹².

Theorem 2. If the assumption (H) is valid, then there exists the positive constants C₃ andC₄ such that

|p^u₂(t)−p^v₂(t)| ≤C3ku−vk₂, ∀t∈[0, T], kI⁰(u)−I⁰(v)k₂ ≤C₄ku−vk₂,

for any u, v∈M_I(u0). Proof. (i) First, from

Lxk(x^u₁, x^u₂, u, t)−Lxk(x^v₁, x^v₂, v, t) =

=−2u−f(x^u₁, x^u₂, t)∂f

∂x_k(x^u₁, x^u₂, t) + 2v−f(x^v₁, x^v₂, t) ∂f

∂x_k(x^v₁, x^v₂, t), using the assumption (H) we deduce that

L_xk(x^u₁, x^u₂, u, t)−L_xk(x^v₁, x^v₂, v, t)≤

≤2C|u(t)−v(t)|+ 4C |x^u₁(t)−x^v₁(t)|+|x^u₂(t)−x^v₂(t)|

≤2C|u(t)−v(t)|+ 4C(C1+C2)ku−vk₂. Then, from

p^u₂(t)−p^v₂(t) =

= Z t

0

n

Lx2(x^u₁(s), x^u₂(s), u(s), s)−Lx2(x^v₁(s), x^v₂(s), v(s), s)−

−(t−s)L_x1(x^u₁(s), x^u₁(s), u(s), s)−L_x1(x^v₁(s), x^v₁(s), v(s), s)^ods−

− t T

Z T 0

n

L_x2(x^u₁(s), x^u₂(s), u(s), s)−L_x2(x^v₁(s), x^v₂(s), v(s), s)−

−(T −s)Lx1(x^u₁(s), x^u₁(s), u(s), s)−Lx1(x^v₁(s), x^v₁(s), v(s), s)^ods,

using the above inequalities we have

|p^u₂(t)−p^v₂(t)| ≤

≤2 Z T

0

2C|u(s)−v(s)|+ 4C(C1+C2)ku−vk₂ds+

+²

√ 3 3 T³²

Z _T

0

2C|u(s)−v(s)|+ 4C(C₁+C₂)ku−vk₂²ds ¹₂

≤C₃ku−vk₂, withC₃ = 4C√

T+ 8C(C₁+C₂) + ⁴

√ 3

3 T³²C^p2 + 8(C₁+C₂).

(ii) From the equality I⁰(u)(t)−I⁰(v)(t) =

=Lu(x^u₁, x^u₂, u, t)−p^u₂(t)−Lu(x^v₁, x^v₂, v, t)−p^v₂(t)

= 2u(t)−v(t)−2f(x^u₁, x^u₂, t)−f(x^v₁, x^v₂, t)−p^u₂(t)−p^v₂(t) we obtain

I⁰(u)(t)−I⁰(v)(t)≤2|u(t)−v(t)|+ (2C+C₃)ku−vk₂. Hence

kI⁰(u)−I⁰(v)k= Z T

0

I⁰(u)(t)−I⁰(v)(t)

2dt

1 2

≤ Z _T

0

h2u(t)−v(t)+ (2C+C3)ku−vk₂ⁱ²dt ¹₂

≤C4ku−vk₂

and C4= 2 + (2C+C3)T.

LetU be a Hilbert space andJ :U →Ra Gˆateaux differentiable functional.

We shall establish an adequate convergence theorem for the gradient method used to solve the optimization problem

minu∈UJ(u).

Theorem 3. Let u⁰ ∈U. If

(1) J is Gˆateaux differentiable and bounded below;

(2) There existsL >0 such that

kJ⁰(u)−J⁰(v)k ≤Lku−vk, for anyu, v∈M_J(u⁰₎=u∈U :J(u)≤J(u⁰) ;

then there exists δ∈(0,_L¹) such that the sequence (u^k)_k∈N, defined by

u^k+1 =u^k−µ_kJ⁰(u^k), with µ_k∈E_δ = [δ,_L² −δ], has the properties:

a) The sequence (J(u^k))_k∈N is convergent;

b) lim_k→∞J⁰(u^k) = 0.

Proof. First we prove that there exists δ ∈ (0,_L¹) such that for any u ∈ M_J(u0), for anyµ∈E_δ and for any t∈[0, µ] we have u−tJ⁰(u)∈M_J(u).

Let us suppose, by contrary, that for any δ ∈ (0,_L¹) there exists u1 ∈ M_J(u⁰₎, µ₁ ∈E_δ and t₁∈[0, µ₁] such that

u₁−t₁J⁰(u₁)6∈M_J(u1) ⇔J(u₁−t₁J⁰(u₁))> J(u₁).

Obviously J⁰(u₁)6= 0. From

λ→0lim

1 λ

J(u₁−λJ⁰(u₁))−J(u₁)=−kJ⁰(u₁)k² <0

it follows that there exists µ₂ such that for any µ∈[0, µ₂]J(u₁−µJ⁰(u₁))<

J(u1). Necessarilyµ2 < t1. The continuity of the functiont7→J(u1−tJ⁰(u1)) implies that there exists t₂ ∈[µ₂, t₁] such that J(u₁−t₂J⁰(u₁)) =J(u₁) and for anyt∈[0, t2], u1−tJ⁰(u1)∈M_J(u1).

The following relations are then valid

0 =J(u₁−t₂J⁰(u₁))−J(u₁)≤ ^Lt₂²² −t₂kJ⁰(u₁)k<0, which are contradictory.

Consequently, the assertions of the theorem follows from the inequalities J(u^k+1)−J(u^k)≤ ^Lµ₂^2k −µ_kkJ⁰(u^k)k ≤ ^Lδ₂² −δkJ⁰(u^k)k, ∀k∈N.

Because the functionalI is Gˆateaux differentiable, bounded below and satis- fies the Lipschitz property (Theorem 3.2), as a consequence of the Theorem 3.3 we obtain the following result.

Theorem 4. If the hypothesis (H) is valid then the sequence (u^k)_k∈N con- structed by the gradient method (17) to solve the NTPBVP (1)–(3) has the following properties:

(1) The sequence (I(u^k))_k∈N is convergent;

(2) lim_k→∞I⁰(u^k) = 0.

4. IMPLEMENTATION OF THE METHOD

Let n ∈ N^∗. On [0,1] we consider the mesh 0 = t₀ < t₁ < . . . < t_n = 1 whereti =ih, i= 0,1, . . . , nand h= 1/n. Let

u^k_h = (u^k₀, u^k₁, . . . , u^k_n), x^k_1,h = (x^k_1,0, x^k_1,1, . . . , x^k_1,n), x^k_2,h = (x^k_2,0, x^k_2,1, . . . , x^k_2,n), p^k_2,h = (p^k_2,0, p^k_2,1, . . . , p^k_2,n)

be the discretization of the functionsu^k, x^u₁^k, x^u₂^k and p^u₂^k respectively, at the pointsti, i= 0,1, . . . , n.

Using the formulas (18), (19), (20) and (6) x^k_1,h, x^k_2,h, p^k_2,h and I(u^k_h) were computed with the trapezoidal rule of integration.

Ifs^k_h = (s^k₀, s^k₁, . . . , s^k_n) are defined by

−s^k_i = 2u^k_i −f(x^k_1,i, x^k_2,i, ti)−p^k_2,i i= 0,1, . . . , n,

then using an algorithm of one dimensional optimization based on a parabolic interpolation, it is findµk as

I(u^k_h+µ_ks^k_h) = minI(u^k_h+µs^k_h) :µ≥0 . The next approximation is given by

u^k+1_i =u^k_i +µks^k_i i= 0,1, . . . , n.

The stopping condition is given by Xⁿ

i=0

(u^k+1_i −u^k_i)^21/2 < = 0.001.

5. NUMERICAL EXAMPLES

Example 1. Consider the equation

x¨= expx, x(0) =x(1) = 0 with the solution

x(t) = ln 2 + 2 ln c/cos^c(t−0.5)₂ , wherec≈1.3360656. In this casef(x₁, x₂, t) =e^x1.

The results are presented in Table 1. On the other hand, the value of the cost functional I(u^k_h) and the error

e_k=

n

X

i=0

[x^k_1,i−x(t_i)]^21/2

are presented in Table 2.

Example 2. Consider the equation

−d dt

h(x²+ 0.1)dx dt

i+x= 10t⁴−20t³+ 11t²+ 0.2, x(0) =x(1) = 0 with the solution x(t) = t−t² (Sokolowski J., Matsumura T., Sakawa Y., [12]). In this case

f(x₁, x₂, t) = x1−2x1x²₂−10t⁴+ 20t³−11t²+t−0.2

x²₁+ 0.1 .

The results are presented in Table 3 and Table 4, respectively.

Remark. The discretization was done with n = 10(h = 0.1). The initial approximations were taken u⁰_i = 0,i= 0,1, . . . , n.

Table 1. The discrete solution. Table 2. The evolution of the cost functional.

tj x^k_1,j x1(tj) |x^k_1,j−x1(tj)|

0 0 0 0

.1 −.0414 −.0414 .2750E−4 .2 −.0732 −.0732 .5038E−4 .3 −.0957 −.0958 .6596E−4 .4 −.1092 −.1092 .7502E−4 .5 −.1136 −.1137 .7799E−4 .6 −.1092 −.1092 .7502E−4 .7 −.0957 −.0958 .6596E−4 .8 −.0732 −.0732 .5038E−4 .9 −.0414 −.0414 .2750E−4

1 0 0 0

k I(u^k) ek

1 1.0000000 .26327643E+0 2 .46221580E−2 .79034139E−2 3 .19229175E−4 .13256727E−2 4 .89363358E−7 .14495543E−3 5 .41923597E−9 .18614450E−3

Table 3. The discrete solution. Table 4. The evolution of the cost functional.

tj x^k_1,j x1(tj) |x^k1,j−x1(tj)|

0 .0000 .0000 0

.1 .0899 .0900 .9783E−5 .2 .1600 .1600 .1367E−4 .3 .2100 .2100 .1578E−4 .4 .2400 .2400 .1712E−4 .5 .2500 .2500 .1761E−4 .6 .2400 .2400 .1712E−4 .7 .2100 .2100 .1578E−4 .8 .1600 .1600 .1367E−4 .9 .0899 .0900 .9783E−5

1 .0000 .0000 0

k I(u^k) ek

1 .15649330E+2 .57732140E+0 2 .67351863E+0 .23050035E−1 3 .22708556E+0 .68533936E−1 4 .84461831E−1 .84342953E−2 5 .30519258E−1 .26329220E−1 6 .11536148E−1 .33611155E−2 7 .43714834E−2 .10240440E−1 8 .16967369E−2 .14039690E−2 9 .66307694E−3 .40626219E−2 10 .26265512E−3 .58864895E−3 11 .10459776E−3 .16319149E−2 12 .41961555E−4 .24410568E−3 13 .16892649E−4 .65968634E−3 14 .68264829E−5 .10047550E−3 15 .27639249E−5 .26765030E−3 16 .11212327E−5 .41146389E−4 17 .45530150E−6 .10879322E−3 18 .18506415E−6 .16800688E−4

REFERENCES

[1] Ascher, U., Christiansen, J. and Russel, R. D.,COLSYS-A collocation code for boundary value problems, Proceedings of Working Conference for Codes for Boundary Value Problems in ODE’s, Houston, Texas, 1978.

[2] Ascher, U., Christiansen, J. and Russel, R. D., A collocation solver for mixed order systems of boundary value problems, Math. Comp.,33, pp. 659–679, 1979.

[3] Fedorenko, R. P., Approximate Solutions for Optimal Control Problem, Nauka, Moskva, 1973 (in Russian).

[4] Goh, C. J. and Teo, K. L., Control parametrization: a unified approach to optimal control problems with general constraints, Automatica,24, no. 1, pp. 3–18, 1988.

[5] Goh, C. J.andTeo, K. L.,MISER: a FORTRAN program for solving optimal control problems, Adv. Eng. Software,10, no. 2, pp. 90–99, 1988.

[6] Keller, H. B., Numerical Solution of Two Point Boundary Value Problems, SIAM Regional Conf., Ser. Appl. Math.,24, SIAM, Philadelphia, 1976.

[7] Klessig, R.andPolak, E.,An adaptive precision gradient method for optimal control, SIAM J. Control, 11, no. 1, pp. 80–93, 1973.

[8] Marzulli, P., Global error estimates for the standard parallel shooting method, J.

Comput. Appl. Math.,34, pp. 233–241, 1991.

[9] Miele, A., Recent advances in gradient algorithms for optimal control problems, J.

Optim. Theory Appl.,17, pp. 361–430, 1975.

[10] Polak, E.,Computational Methods in Optimization, Academic Press, New York, 1971.

[11] Scheiber, E.,Numerical solution of a nonlinear two point boundary value problem by optimal control methods, Bull. Univ. Bra¸sov,30, ser. C, pp. 51–56, 1988.

[12] Sokolowski, J., Matsumura, T.andSakawa, Y.,Numerical solution of a nonlinear two point boundary value problem by an optimization technique, Control and Cybernet- ics,11, nos. 1–2, pp. 41–56, 1982.

[13] Teo, K. L., Goh, C. J. and Wong, K. H., A Unified Computational Approach to Optimal Control Problems, Longman Scientific & Technical, New-York, 1991.

Received by the editors: November 25, 1999.