Solution Methods of Optimal Control Problems
Robert Baier
, Christof Buskens y
, Ilyes Assa Chahma z
, Matthias Gerdts x
13th Deember 2004
Abstrat
A numerial method for the approximation of reahable sets of linear
ontrol systems is disussed. The method is based on the formulation of
suitableoptimal ontrolproblems with varying objetive funtion, whose
disretizationbyRunge-Kutta methods lead to nitedimensional onvex
optimizationproblems. Itturnsoutthattheorderofapproximationforthe
reahablesetdependsonthepartiularhoieoftheRunge-Kuttamethod
inombinationwiththeseletionstrategyusedforontrolapproximation.
For an inappropriate ombination the expeted order of onvergene an
not be ahieved in general. The method is illustrated by two examples
usingdierent Runge-Kutta methods and seletion strategies and allows
to estimatetheorder of onvergenenumerially.
Keywords: optimalontrol, approximationof reahable sets,diret solution
methods, order of onvergene
Department of Mathematis, University of Bayreuth, 95440 Bayreuth, Germany,
Robert.Baieruni-bayreuth.de
y
DepartmentofMathematis,UniversityofBremen,28344Bremen,Germany
z
Cetelem Bank GmbH, Shwanthalerstrasse 31, 80336 Munhen, Germany,
a.hahmaetelembank.de
x
Department of Mathematis, University of Bayreuth, 95440 Bayreuth, Germany,
Matthias.Gerdtsuni-bayreuth.de
Thesubjetofthispaperisthedesriptionofanalgorithmfortheapproximation
of reahable sets of linearontrolproblems. The problemof determining onvex
reahable sets an be equivalently desribed by innitely many optimal ontrol
problems, where the objetive funtion is adapted. By hoosing only nitely
many diretionsapproximationsof reahable sets anbeobtained. Theouring
optimalontrolproblems are not solved theoretially by use of the Pontryagin's
maximumprinipleasin[38℄butnumeriallybysuitabledisretizationmethods.
This allows to treat also time dependent linear problems and even nonlinear
ones. Non-polyhedral ontrol regions an be treated as nonlinear inequalities
andequalities. Results onerningthe onvergene of disretizedoptimalontrol
problems an be found in [30℄, [10℄ and the referenes stated therein.
Inthis ontext, thepartiular hoieof theseletionstrategyused forontrol
approximationturnsouttoberuialfortheorderofonvergeneanddependson
thehoieoftheRunge-Kuttashemeusedforthedisretizationoftheunderlying
dierentialequations. Inordertoillustratethis dependenyseveralRunge-Kutta
methods with dierent seletionstrategies (pieewise onstant, pieewise linear,
independent seletion) are disussedin moredetail fortwo illustrativeexamples.
By this approah umbersome set operations (like Minkowski sums, unions
of sets, ...) an be avoided and lead to known optimization methods, whih
in addition yield not only the endpoints of optimal trajetories, but the entire
trajetory inluding the orresponding optimal ontrol. Furthermore, this ap-
proah isuseful forlinear ontrolproblems with ontrolregionsformulatedwith
nonlinearrestritions(see (7))and innonlinearontrolproblemsyielding onvex
reahable sets, too. However, the lose onnetion between set-valued analysis
andoptimalontrolisshowninSetion3. Aomparisonwithset-valuedmethods
asin [12, 4,3, 41, 8℄ isbeyond the sope of this paper.
Methodsforlineardierentialinlusionsbasedonset-valuedquadraturemeth-
ods or set-valued Runge-Kutta methods are mentioned in [3℄ as well as other
methods, e.g. estimation methods for reahable sets (f. [15℄) and ellipsoidal
methods(f. [23℄ foranoverview). Newerdevelopmentsofthesemethodsahieve
inner approximations ([24℄, [26℄) and outer approximations [25℄ of the reahable
set (see also [4℄).
The problem of the approximation of reahable sets appears inseveral disi-
plines: ontrol theory, ordinary dierential equations with unertainties or with
disontinuities in the state, neessary onditions for a minimum in nonsmooth
analysis, dierential games and viability theory, f. [5℄, [1℄, [33℄, [14℄. The on-
vexity ofthese reahablesets an beguaranteed for lineardierentialinlusions,
but may alsoappear for nonlinear problems.
The paper is organized as follows. In Setion 2 basi notations and proper-
ties of reahable sets are summarized. Basi fats on the desription of onvex
sets and arithmeti set operations are introdued and form the basis for the re-
dened, whih are used to measure the speed of onvergene w.r.t. the opti-
mal value and the optimal trajetory, respetively. In Setion 3 the problem of
alulating the boundary of the reahable set is reformulated as innitely many
optimalontrolproblems whih dieronlyin theobjetivefuntion. Theseopti-
malontrolproblemsaredisretizedbyuse ofexpliitRunge-Kuttamethodsand
suitableontrolapproximationsresultinginnitedimensional (linear/nonlinear)
optimizationproblems. Herein,several approximationlassesfortheontrollead
to dierent seletion strategies in the disretization. The setion ends with a
formulationof the proposed methodfor the approximation ofreahable sets and
its implementation. Several ombinations of Runge-Kutta methods and sele-
tion strategies are disussed in Setion4 with illustrativeexamples. Tables with
onvergene resultsand visualizations ofreahable sets are inluded. Finally, an
outlinefor further researh onludes the paper.
2 Notation
In this setion,some introdutory denitions and results are olleted.
The basi underlyingproblem isthe following ontrol problem:
Problem 2.1 LetA():R n
!R nn
andB():R m
!R m n
betwoL
1
-integrable
matrix funtions.
LetU R m
beanonempty,onvexompatsetandI :=[t
0
;t
f
℄bea realinterval.
For a given ontrol funtion u : I ! R m
with u() 2 L
1 (I;R
m
) we are looking
for a solution x()2W 1;1
(I;R n
) of the dierential equation
_
x(t)=A(t)x(t)+B(t)u(t) (a.e. t2I); (1a)
x(t
0 )=x
0
; (1b)
u(t)2U (a.e. t2I). (1)
Denition 2.2 Let usstudy Problem 2.1 and let t2I. Then,
R(t;t
0
;x
0
):=fy 2R n
j9u() ontrol funtion and 9x() orresponding
solution of Problem 2.1 with x(t)=yg
is alled the reahable set of the orresponding ontrolproblem for the time t.
In 1965, Aumann disovered the onvexity of the set-valued integral in [2℄
whiheasilyleadstotheonvexityofthereahablesetforlinearontrolproblems.
Proposition 2.3 In Problem 2.1, the reahable set R(t;t
0
;x
0
) is onvex, om-
pat and nonempty for every t2I.
SomenotationsfromConvexAnalysisarerealledwhihareneessary forthe
explanation of the algorithmdesribed later.
Denition 2.4 Denote by C(R n
) the setof allnonempty onvex ompat sets in
R n
and let C 2C(R n
) and l2R n
.
Then,
Æ
(l;C):=max
2C l
>
isthe support funtion of C in diretion l and
Y(l;C):=f2Cjl
>
=Æ
(l;C)g
is the set of supporting points of C in diretion l.
We need the following property of support funtions:
Lemma 2.5 Let C = C
1
C
2
2 C(R n
) with onvex sets C
i
R
n
i
, n
i 2
f1;:::;ng, i = 1;2, and n
1 +n
2
= n. Then, for given l = (l
>
1
;l
>
2 )
>
2 R n
with l
i 2R
n
i
, i=1;2, we have:
Æ
(l;C)=Æ
(l
1
;C
1 )+Æ
(l
2
;C
2 ):
Proof: see e.g. [19, xV, Disussionafter Remark 3.3.6℄
Support funtions resp. supporting points desribe fully a onvex ompat
set.
Proposition 2.6 Let C 2C(R n
). Then,
C =
\
kl k2=1
fx2R n
jl
>
xÆ
(l;C)g; C = [
kl k2=1
Y(l;C);
C =o(
[
kl k
2
=1
fy(l;C)g) with arbitrary y(l;C)2Y(l;C);
where C denotes the boundary of C and o () denotes the onvex hull of a set.
Proof: see e.g. [19, xV.,Theorem 2.2.2℄and [19, xV.,Proposition 3.1.5℄.
The last equation follows easily, if one estimates the support funtion of the
right-hand side in diretion by
>
y(;C)=Æ
(;C) frombelow.
A ommon arithmeti operationson sets is the salar multipliation and the
Minkowski sum whihare realledhere.
Denition 2.7 Let C;D2C(R ), 2R and A2R . Then,
C :=fj2Cg
denes the salar multipliation,
AC :=fAj2Cg
the image of C under the linear map x7!Ax and
C+D:=f+dj2C;d2Dg
the Minkowski sum.
We needthe followingtheoretialresultwhihstates onvexity and ompat-
ness of the set operationsdened above.
Lemma 2.8 Let C;D 2 C(R n
), 2 R and A 2 R mn
. Then, C and C +D
are elements of C(R n
) and AC is an element of C(R m
). Furthermore,
Æ
(l;C)=Æ
(l;C); Y(l;C)=Y(l;C) (if 0);
Æ
(;AC)=Æ
(A
>
;C); Y(;AC)=AY(A
>
;C);
Æ
(l;C+D)=Æ
(l;C)+Æ
(l;D); Y(l;C+D)=Y(l;C)+Y(l;D)
for all l2R n
, 2R m
.
Proof: Toguarantee thatthe operationsgiveresultsinC(R n
)and theequations
on the support funtions see [19, xV, Theorem 3.3.3(i) and Proposition 3.3.4℄.
The equations on the supporting set follow immediately from alulus rules on
the subdierential in [19, xVI, Theorem 4.1.1 and equation (3.1)℄ and [32, The-
orem 23.9℄, sine [19, xVI, Proposition 2.1.5 and equation (3.1)℄ onnets the
subdierentialof the support funtion and the supporting set.
Denition 2.9 Let C;D2C(R n
). Then,
d(C;D):=max
2C min
d2D
k dk
2
;
d
H
(C;D):=maxfd(C;D);d(D;C)g
are dening the one-sided Hausdordistane resp. the Hausdordistane of the
two sets.
The Demyanov distane between two sets is dened as
d
D
(C;D):= sup
l 2T
C
\T
D
ky(l;C) y(l;D)k
2
;
whereT
C
isdenedas setof allnormed diretionsin R n
forwhih thesupporting
set Y(l;C) onsists of only one point y(l;C) (T
D
is dened analogously for the
set D). T
C
and T
D
are subsets of the unit sphere of full measure.
following result forthe Hausdor-distane:
Proposition 2.10 Let C;D2C(R n
). Then,
d
H
(C;D)= max
kl k2=1 jÆ
(l;C) Æ
(l;D)jd
D
(C;D):
Proof: see e.g. [19, xV, Theorem 3.3.8℄and [9,Lemma 4.1℄
3 New Method for the Approximation of Reah-
able Sets
3.1 Computation of the Reahable Set by Optimal Con-
trol
Sine we know from Proposition 2.3 that the reahable set for problem 2.1 is
onvex, it issuÆient toalulate merely the boundary of the reahable set.
Proposition 2.6 gives a motivation to alulate at least one support point
(whihliesautomatiallyattheboundary)ofthereahablesetindiretionl 2IR n
with klk
2
= 1. Note that even in the ase that the reahable set is not stritly
onvex and the set of supportingpointsisa (n 1)-dimensionalfae,for a xed
diretionl, one supporting point inthis diretionis suÆient to reonstrut the
reahable set.
Thus, to alulate a supporting point x(t
f
) on the boundary of the reah-
able set R(t
f
;t
0
;x
0
) in a xed diretion l we have to nd an admissible ontrol
funtion u(t) 2 U that maximizes the funtional y 7! l
>
y (resulting in the sup-
portfuntionÆ
(l;R(t
f
;t
0
;x
0
))asoptimal value). This onstitutes the following
speial optimal ontrol problemof Mayer type:
(OCP
l )
8
<
:
Maximize l
>
x(t
f )
w.r.t. u2L 1
([t
0
;t
f
℄;IR m
);x2W 1;1
([t
0
;t
f
℄;IR n
)
x() orrespondingsolution tou() for(1a){(1).
We denote the optimal solution of (OCP
l
) by x
?
(t;l) and u
?
(t;l), where the
argument l indiates the dependeny of the diretionl.
AsalreadymentionedinProposition2.6,theonvexityandompatnessofthe
reahablesetguaranteedbyProposition2.3leadstotheequivalentrepresentation
by onsidering supporting points inall diretions l 2IR n
,klk
2
=1:
R(t
f
;t
0
;x
0
)=ofx
?
(t
f
;l)jl 2IR n
;klk
2
=1g:
timal Control Problems
In general, for omplex problems neither we an ompute a solution of (OCP
l )
analytially nor for all diretions l. Hene, we suggest to approximate (OCP
l )
numeriallyandonsideronlyanitenumberofdiretionsl
i
,i=1;:::;M :=N
l .
This yieldsanapproximation
R
M (t
f
;t
0
;x
0
)R(t
f
;t
0
;x
0 )
of the reahable set whih willbespeied hereafter.
Forthe momentlet l be xed with klk
2
=1.
ForN
t
2IN;N
t
2 weintroduea grid with grid points
t
i
=t
0
+ih2[t
0
;t
f
℄;i=0;1;:::;N :=N
t
;h= t
f t
0
N
t
: (2)
The ontrol funtion u(t) is disretized on eah subinterval [t
i
;t
i+1
℄ by the ap-
proximation
u (i)
app
(t;^u); t2[t
i
;t
i+1
℄;
where ^u = (u
0
;u
1
;:::;u
P 1 )
>
2 U P
is a nite dimensional vetor parametriz-
ing the seletion strategy for the ontrol in the following expliit Runge-Kutta
sheme.
Letusrst deneexpliitRunge-Kuttashemesbeforewewilldisussparti-
ularstrategiesforthe approximationofthe ontrolinmoredetails. Eahexpliit
Runge-Kuttasheme an beharaterized by its Buther array:
1
0 0
2
21
0 0
.
.
. .
.
. .
.
. .
.
. .
.
.
s
s1
s;s 1 0
1
s 1
s
For a given ontrol approximation u (i)
app
(t;u)^ on [t
i
;t
i+1
℄ a state approximation
x
app
(t;^u) isobtained viaan expliit s-step Runge-Kutta disretizationsheme:
x
app (t
i+1
;u)^ = x
app (t
i
;u)^ +h(x
app (t
i
;u);^ ^u;h); i=0;1;:::;N
t 1;
x
app (t
0
;u)^ = x
0
(3)
and
(x
app (t
i
;^u);^u;h):=
s
X
j=1
j
A(t
i +
j h)
(j)
i+1 +B(t
i +
j h)u
(i)
app (t
i +
j h;u)^
;
(j)
i+1 :=x
app (t
i
;u)^ +h j 1
X
k=1
jk
A(t
i +
k h)
(k)
i+1 +B(t
i +
k h)u
(i)
app (t
i +
k h;u)^
:
jk j j
[7℄.
Letus now onsider examples forseletion strategies used inSetion 4.
(i) Continuous and pieewise linearapproximation:
u (i)
app
(t;u )^ :=u
i +
t t
i
h (u
i+1 u
i
) for t2[t
i
;t
i+1
℄;i=0;1;:::;N 1;
withP =N +1.
(ii) Pieewise onstant approximation:
u (i)
app
(t;^u):=u
i
for t2[t
i
;t
i+1
℄;i=0;1;:::;N 1;
withP =N.
(iii) Independent seletions atintermediate grid points t
i +
j h:
u (i)
app (t
i +
j
h;^u):=u
is+j 1
; i=0;1;:::;N 1; j =1;:::;s; (4)
with P =sN. Here, eah grid pointreates a new independent seletion
foreahsubinterval. FormodiedEuler'smethod(seeSetion4andFigure
4 in Example 4.2)
1
= 0,
2
= 1
2
so that two independent seletions u
2i
and u
2i+1
are hosen fromU for this method ineah subinterval[t
i
;t
i+1
℄.
For Heun's method (see Setion 4 and Figure 3 in Example 4.2)
1
= 0,
2
=1sothattwoindependentseletionsu
2i andu
2i+1
arealsohosenfrom
U forthismethodineahsubinterval[t
i
;t
i+1
℄,althought
i +
2 h=t
i+1 +
1 h
fori=0;:::;N 1.
Please notie, that further seletionstrategies are possible, e.g. independent se-
letions with additional ontinuity onstraints at the inner grid points t
i , i =
1;:::;N 1,oradditionalequalityonstraints atthose intermediategrid points
t
i +
j
h where dierentindies j produe the same intermediate grid point(i.e.,
pointswhere
j
=
k
with j 6=k).
Thus, by this disretization the innitedimensional optimal ontrol problem
(OCP
l
) is approximated by the nite dimensionalonvex programming problem
(CP 1
l )
8
>
>
>
>
>
>
<
>
>
>
>
>
>
:
Maximize l
>
x
app (t
Nt
;^u)
w.r.t. u^ 2U P
subjet to
x
app (t
i+1
;^u) = x
app (t
i
;^u)+h(x
app (t
i
;u);^ u;^ h);
i=0;1;:::;N
t 1
x
app (t
0
;^u) = x
0
;
^
u 2 U
P
: (?)
Notie,that^uimpliitlydenes aontrolapproximationu (i)
app
(;^u)oneahsubin-
terval [t
i
;t
i+1
℄,ompare the examples(i)-(iii).
We denote the optimalsolution of (CP
l
) by u^ .
Iftheonditions(?) anbewrittenwithanitenumberofaÆneinequalities,
(CP 1
l
) isa linearprogramming problemand alled(LP 1
l
),otherwise anonlinear
(onvex) programmingproblem.
In the sequel, we investigate the simplestase, the Euler's methodIn the se-
quel,weinvestigatethesimplestase,theEuler'smethodwithpieewiseonstant
ontrolapproximation,sine it isthen easier possible toderive expliit solutions
for the nite dimensional problems (CP 1
l
). Nevertheless, every expliit Runge-
Kutta methods with the seletion strategies (i){(iii) will give a similar (more
ompliated) representation. The expliit formulae for the solution stress the
strong onnetionto set-valued methods e.g. in [12, 4, 41℄via supportfuntions
resp. supporting points.
In the ase of Euler, (3)redues to
(x
app (t
i
;^u);^u;h)=A(t
i )x
app (t
i
;^u)+B(t
i )u
i :
The reursive evaluationin (3)for Euler's method yields
x
app (t
Nt
;u)^ = Nt 1
Y
i=0 Q
i
!
x
0 +h
Nt 1
X
k=0
Nt 1
Y
i=k+1 Q
i
!
B
k u
k
(5)
withQ
i
:=I+hA(t
i ), B
k
:=B(t
k
) and the nn-identitymatrix I. The matrix
produt Q
is dened as
j
Y
i=k Q
i :=Q
j Q
j 1 Q
k :
Introduing this expression for x
app (t
f
;^u) in(LP 1
l
) yieldsthe linear program
(LP 2
l )
8
>
<
>
:
Maximize l
>
N
t 1
X
k=0 N
t 1
Y
i=k+1 Q
i
!
B
k u
k
!
subjet to u
k
2U; k =0;1;:::;N
t 1:
Note that this linear program has the same solution u^ as (LP 1
l
), whereas the
optimalobjetivefuntionvaluesare dierent,sinewenegletedonstantterms.
To ompute the objetive funtion in(LP 2
l
)very eÆiently we introdue ad-
ditionalartiialvariables
>
Nt
:= l
>
;
>
i
:=
>
i+1 Q
i
=
>
i+1 +h
>
i+1 A
i :
These artiialvariables are alulated bakward in time and orrespond to the
disretized adjointvariable of the optimalontrol problem(OCP
l ).
Then, (LP
l
) an be replaed by
(LP 3
l )
8
>
<
>
:
Maximize Nt 1
X
k=0
>
k+1 B
k u
k
subjet to u
k
2U; k =0;1;:::;N
t 1:
Lemma 2.8gives us
Nt 1
X
k=0 Æ
(
k+1
;B
k U)=
Nt 1
X
k=0 Æ
(B
>
k
k+1
;U)
as optimal value of (LP 3
l
) and hene, (u
0
;u
1
;:::;u
N
t 1
) with the supporting
pointsu
k
2Y(B
>
k
k+1
;U) as one solution.
In the speial of box onstraints, that is U = fu 2 IR m
j u u ug, we
deneS
>
k
:=(S 1
k
;:::;S m
k
):=
>
k+1 B
k 2IR
m
. Sine the objetive funtion
Nt 1
X
k=0 S
k u
k
= Nt 1
X
k=0 m
X
j=1 S
j
k u
j
k
ismaximized,if eahtermS j
k u
j
k
ismaximized,the solutionof (LP 3
l
)isgiven by
u j
k
= 8
>
>
>
<
>
>
>
: u
j
; if S
j
k
<0;
u j
; if S
j
k
>0;
arbitrary; else:
forj =1;:::;m; k=0;:::;N
t 1.
3.3 Disrete reahable sets
Disrete reahable sets are the reahable sets of the disretized equations and
ouldbedened asendpoints ofdisrete solutionsof the following problem.
Given thedata inProblem2.1,the disretizedproblemdependsonthe hoie
ofthesetU
app
ofalldisretizedontrolfuntionsandontheRunge-Kuttasheme.
Problem 3.1 For a time disretization (2) with step-size h= t
f t
0
Nt
and a given
disretized ontrol funtion u
app
(;u)^ we are looking for a solution x
app
(;^u) at
the grid-points t
i
, i=0;1;:::;N
t , with
x
app (t
i+1
;^u)=x
app (t
i
;^u)+h(x
app (t
i
;^u);u;^ h) (6a)
for i=0;1;:::;N
t 1;
x
app (t
0
;^u)=x
0
; (6b)
u
i
2U; i=0;1;:::;N
t
; (6)
u
app
(;^u)2U
app :
f0;1;:::;N
t
g. Then,
R
N (t
i
;t
0
;x
0
):=fy2R n
j9u
app
(;u)^ disretized ontrol funtion and
9x
app
(;u)^ orresponding solution of Problem 3.1
with x
app (t
i
;^u)=yg
isalledthe disretereahablesetoftheorrespondingdisretizedontrolproblem
for the time t
i .
Thedenitionaboveshowsthateahoptimizerofproblem(CP 1
l
)(resp.therefor-
mulation(LP 3
l
)) isasupporting pointof thedisretereahableset R
N (t
f
;t
0
;x
0 )
in diretion l. The optimal value of problem (CP 1
l
) oinides with the support
funtionÆ
(l;R
N (t
f
;t
0
;x
0
)). Proposition 2.6shows that
R
N (t
f
;t
0
;x
0 )=
\
kl k
2
=1
fx2R n
jl
>
xl
>
x
app (t
f
;u^
?
)g;
R
N (t
f
;t
0
;x
0
)=o(
[
kl k2=1 fx
app (t
f
;^u
?
)g):
In pratie, onlya nite number of dierent normeddiretions l i
, i=1;:::;M,
are hosen.
Proposition 3.3 ConsiderProblem3.1withatimedisretization (2)andleti2
f0;1;:::;N
t
g. Then, the orresponding disrete reahable set is onvex, ompat
and nonempty.
Proof: For a hosen disretized ontrolfuntion u
app
(;^u), the disrete solution
isdened by (5). The disretereahable set oinides with the union of allsuh
disretesolutionsforallfeasibledisretizedontrolfuntions. IntheaseofEuler
and linear approximation of the ontrols, this orresponds to the union over all
vetors u^2R m (N+1)
. Denition 2.7shows that the disrete reahable set
R
N (t
f
;t
0
;x
0 )=
N
t 1
Y
i=0 Q
i
!
x
0 +h
N
t 1
X
k=0 (
N
t 1
Y
i=k+1 Q
i
!
B
k )U
is a saled Minkowski sum of linearly transformed onvex sets U. Lemma 2.8
proves the wanted properties of the disrete reahable set.
3.4 Implementation
In the sequel, we briey disusssome numerialmethods, whih are suitablefor
solving the disretized optimal ontrol problem (CP 1
l
). Of ourse, the hoie
region U. Hene, we restrit the disussion to onvex ontrol regions U dened
by
U =fu2X jg
i
(u)0; i=1;:::;rg; (7)
where X := fu 2 IR m
j e
Au = b; u 0g with a matrix e
A 2 R pm
and the
funtionsg
i
(),i=1;:::;r, ould be eitherlinear or nonlinear.
Remark 3.4 In the ase, that the support funtion or the supporting points of
theonvexontrolsetU areknown,generalontrolregionsU anbeapproximated
in another way. Proposition2.6 suggests to use the approximation
U
\
i=1;:::;M
fx2R m
j i
>
xÆ
( i
;U)g
resp.
U o(
[
i=1;:::;M fy(
i
;U)g) with arbitrary y(
i
;U)2Y( i
;U) :
Herein, theM dierentnormeddiretions i
2R m
shouldbe hosen inan appro-
priateway inorder toapproximatetheunit sphere. Onemethod istoparametrize
them byspherialoordinates and useequidistant partitionson the parameter in-
tervals for the angles (see [3, Subsetion 3.1.2℄).
If the funtions g
i
in (7) are aÆne linear, then problem (CP 1
l
) is a linear
optimization problem and an be solved by the well-known simplex method or
some interior point method, f. [42℄, suitable for linear programs. In the speial
ase of an Euler approximation and U dened by box onstraints only, a very
eÆientmethodis desribed inSetion 3.2.
If the funtions g
i
are onvex and smooth, i.e. atleast ontinuously dieren-
tiable,thenthe resultingproblem(CP 1
l
)isaonvexbut nonlinear programming
problem and the sequential quadrati programming (SQP) method is appropri-
ate provided the funtions g
i
are dened for infeasible points, f. [34℄, [35℄, [18℄.
Alternatively, the method of feasible diretions is appliable, espeially, if the
funtionsg
i
are only dened for admissiblepoints, f. [43℄.
Ifthefuntionsg
i
are onvexbut nonsmooth,the bundlemethodrespetively
the bundletrust region method (BT-method)is suitable,f. [28℄, [31℄, [21℄, [22℄,
[36℄. Inaddition,Kelly'suttingplanemethodisalsoappliable,f. [20℄. Notie,
that the BT-method and the utting plane method are losely related, f. [21℄,
[36℄.
In the sequel we refer to the optimal ontrol problem (OCP
1
), the dierential
equation (1a)-(1b), the ontrol onstraint (1), and the ontrol approximations
disussed in(i)-(iii) inSetion 3.2.
The following Runge-Kuttamethods are used forthe numerialomputation
of reahable sets:
0 0
1
Euler's method
0 0 0
1 1 0
1=2 1=2
Heun's method
0 0 0
1=2 1=2 0
0 1
Modied Euler'smethod
For all numerial experiments the number of diretions M in Remark 3.4 is
hosenas1200???. Forsimpliity,the methodswithdierentseletionstrategies
are tested for time-independent two-dimensional problems (in whih one ould
even alulate a theoretial solution for referene purposes). Nevertheless, the
framework presented before is still valid and the methods ould be used also
in more ompliated problems (time dependent and higher dimensional) met in
pratie.
From Denition 2.9 of the Hausdor distane, it is lear that the approxi-
mationof thereahable set orresponds toa uniformonvergene of the optimal
value funtions, whereas the approximation of trajetories orresponds to the
uniformonvergene of the maximizersand the Demyanov distane.
Example 4.1 (see [39, Example in setion 4℄) Letusonsiderthefollowing
example with n=2, m =1, x
0
=(0;0)
>
, I =[0;1℄, U =[0;1℄, and
A(t)=
0 1
0 0
; B(t)=
0
1
:
In Figure 1 approximations to the reahable set R(1;0;x
0
) are shown, in the
leftpitureapproximationswith Euler'smethodwithpieewiseonstantseletions
are shown (rst order of onvergene), in the right one the orresponding ones
for Heun's method with ontinuous and pieewise linear ontrol approximation
(seondorder ofonvergene)aredepited. Inbothasesthesetwiththesolidline
showsthe referene set (alulatedwith theorresponding method forN =1280).
The dashed lines show the approximations for N =10;20;40for Euler's method
on the left piture (please note the halfening of the distane of the upper right
orner of the sets when the number of subintervals is doubled). At the right one,
the dashed lines show the approximations for N = 1;2;4 for Heun's method (a
smaller number of subintervals are hosen so that one ould still see in Figure 1
a dierene of the orresponding approximations). Please notie the more rapid
Euler's method.
-0.2 0 0.2 0.4 0.6 0.8 1 1.2
0 0.1 0.2 0.3 0.4 0.5
-0.2 0 0.2 0.4 0.6 0.8 1 1.2
0 0.1 0.2 0.3 0.4 0.5
Figure 1: First order ontra seond order approximations to the reahable set
(left: Euler's method with error O(h), right: Heun's methodwith error O(h 2
))
As Veliov explains in [39℄, the onvergene of the trajetory ould not be bet-
ter than O(h) in this example. In Figure 2 the rst order approximations to the
ontrolandtothestate omponents(oordinates x
1 andx
2
) areshownforHeun's
method with ontinuous, pieewiselinear seletions. Again, the referene isom-
puted by the method itself with N = 1280 (solid line) and in dashed lines the
approximations forN =10;20;40. As it islearly seen, theorder of onvergene
is only1.
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1 0
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
0 0.2 0.4 0.6 0.8 1 0
0.1 0.2 0.3 0.4 0.5
0 0.2 0.4 0.6 0.8 1
Figure 2: First order approximations to the ontrol (left) and the state ompo-
nents (middle,right)by Heun's method
Here, the ombination method of set-valued iterated trapezoidal rule together
withHeun'smethodintroduedin [3,4℄ withN =1000000servesas the referene
set R
ref (0;x
0
;). By omparing the dierent values based on the optimal value
funtion resp. the maximizers, the order of onvergene is estimated. The angle
'for the diretion l 2R 2
, in whih the maximum in (8) resp. (9) is attained, is
shown in the most right olumn.
Hausdor estim.
N distane order angle
10 0.05000000 NaN 0.00500
20 0.02500000 1.00000 0.00500
40 0.01250000 1.00000 0.00500
80 0.00625000 1.00000 0.00500
160 0.00312500 1.00000 0.00500
320 0.00156250 1.00000 0.00500
640 0.00078125 1.00000 0.00500
Demyanov estim.
N distane order angle
10 0.13702925 NaN 5.55500
20 0.06806368 1.00953 5.55500
40 0.03392323 1.00461 5.51500
80 0.01731662 0.97012 5.53500
160 0.00861479 1.00727 5.53500
320 0.00426388 1.01465 5.53500
640 0.00209303 1.02657 5.62500
Table I: order of onvergene for Euler's method (left table: approximation of
the reahable set, righttable: approximation of the trajetories).
Table I shows the expeted order of onvergene 1 for reahable set and the
trajetories. As remarked above the Hausdor distane is attained at the upper
right orner. This table shows the approximated values
max
i=1;:::;M jÆ
(l
i
;R(1;0;x
0 )) Æ
(l
i
; b
R
ref (0;x
0
;))j (8)
resp.
max
i=1;:::;M kY(l
i
;R(1;0;x
0
)) Y(l
i
; b
R
ref (0;x
0
;))k
2
(9)
at the hosen diretions l
i
, i=1;:::;M, for the two distanes
d
H
(R(1;0;x
0 );R
N
(1;0;x
0
)) resp. d
D
(R(1;0;x
0 );R
N
(1;0;x
0 )):
Hausdor estim.
N distane order angle
10 0.00124700 NaN 3.09500
20 0.00031111 2.00295 3.12000
40 0.00007788 1.99805 6.27500
80 0.00001947 1.99990 3.14000
160 0.00000488 1.99688 6.26000
320 0.00000122 1.99929 3.14500
640 0.00000030 2.00266 6.22500
Demyanov estim.
N distane order angle
10 0.06636590 NaN 5.55500
20 0.03273184 1.01975 5.55500
40 0.01668369 0.97226 2.40000
80 0.00848003 0.97630 5.53500
160 0.00419649 1.01488 5.53500
320 0.00205473 1.03024 5.53500
640 0.00099208 1.05042 5.62500
Table II: order of onvergene for Heun's method (left table: approximation of
the reahable set, righttable: approximation of the trajetories)
For Heun's method with ontinuous, pieewise linear ontrol approximation,
Table II shows order of onvergene 2 for the reahable set and only order 1 for
the trajetories.
with n=2, m =2, x
0
=(0;0)
>
, I =[0;2℄, U =fx2IR 2
jkxk
2
1g, and
A(t)=
0 1
2 3
; B(t)=
1 0
0 1
:
This example introdues the nonlinear onstraint
u 2
1 +u
2
2 1
for the ontrol variableu=(u
1
;u
2 )
>
.
The seond order approximations to the reahable set R(2;0;x
0
) alulated
byHeun'smethodwith pieewiseonstant ontrols resp. with independentontrol
seletion in t
i and t
i+1
(see (4)) are shown in Figure 3.
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
-1.5 -1 -0.5 0 0.5 1 1.5
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
-1.5 -1 -0.5 0 0.5 1 1.5
Figure 3: seond order approximations to the reahable set for Heun's method
with pieewise onstant ontrol approximation (left) resp. independent ontrol
seletion(right).
The set with the solid line shows the referene set (alulated with the orre-
sponding methodfor N =160)and thedashedlines representthe approximations
for N =5;10;20. At the left piture the onvergene order O(h 2
) an be seen by
studying the boundary of the sets near by y=1.
Both seletion strategies seems to onverge with order 2 whih is assured by
Tables IIIand IV.
Nevertheless, Figure 4 showsthat the hoieof the seletionstrategies forthe
ontrol should depend on the Runge-Kutta method. In Figure 4 the pieewise
onstantseletion strategy isompared with the independent ontrol seletions in
t
i and t
i +
h
2
for modied Euler's method (see (4)). The latter seletion strategy
distane
5 0.10328935 NaN 1.37000
10 0.02307167 2.16250 1.53000
20 0.00521186 2.14625 1.57500
40 0.00123195 2.08086 4.73500
80 0.00029922 2.04164 1.60000
160 0.00007372 2.02105 4.74500
distane trajetory
5 0.37223126 NaN 0.90500
10 0.07159599 2.37825 0.88500
20 0.01535558 2.22112 4.02500
40 0.00355544 2.11066 4.02500
80 0.00085565 2.05493 4.02500
160 0.00020992 2.02719 4.02500
TableIII: Orderof Convergene for Heun's methodwith pieewise onstanton-
trolapproximation.
N Hausdor Order angle
distane
5 0.04517018 NaN 1.72000
10 0.00772443 2.54787 4.23500
20 0.00203009 1.92789 4.30000
40 0.00051385 1.98211 4.33500
80 0.00012897 1.99429 1.21000
160 0.00003229 1.99784 1.22000
N Demyanov Order angle
distane trajetory
5 0.16781544 NaN 1.18500
10 0.04611042 1.86371 0.87500
20 0.01077148 2.09788 4.01500
40 0.00257389 2.06520 0.87500
80 0.00062808 2.03492 0.87500
160 0.00015506 2.01802 4.01500
Table IV: Order of Convergene for Heun's method with independent seletion
strategy(iii).
destroys order of onvergene 2 of the Runge-Kutta method. This is veried in
the Tables V (order O(h 2
)) and VI (only order O(h)) for the onvergene to the
reahable set and the trajetories.
N Hausdor Order angle
distane
5 0.10328935 NaN 1.37000
10 0.02307167 2.16250 1.53000
20 0.00521186 2.14625 1.57500
40 0.00123195 2.08086 4.73500
80 0.00029922 2.04164 1.60000
N Demyanov Order angle
distane trajetory
5 0.37223121 NaN 0.90500
10 0.07159599 2.37825 0.88500
20 0.01535559 2.22112 4.02500
40 0.00355571 2.11056 4.02500
80 0.00085566 2.05503 0.88500
Table V: Order of Convergene for the modied Euler's method with pieewise
onstantontrolapproximation.
distane
5 0.83583108 NaN 4.03000
10 0.33319435 1.32685 0.85500
20 0.15333206 1.11970 5.34000
40 0.07575471 1.01725 5.36000
80 0.03762644 1.00959 2.22500
distane trajetory
5 1.03202096 NaN 0.73000
10 0.36562913 1.49702 3.85000
20 0.16060144 1.18690 3.76000
40 0.07933801 1.01740 4.72000
80 0.03952243 1.00534 4.72000
Table VI: Orderof Convergene for the modiedEuler's method with freesele-
tion.
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
-1.5 -1 -0.5 0 0.5 1 1.5
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
Figure 4: approximations to the reahable set for N = 160 (solid) and N =
5;10;20(dashed) omputedby modied Euler'smethodwithpieewise onstant
(left)resp. independent seletionstrategy (iii)(right).
It is known that set valued quadrature methods in [4℄ ould lead to a order of
onvergene greaterthantwo,ifthe problemsatisesadditionalsmoothnesson-
ditions,f. [3℄. In this ase, seletion strategies with pieewise onstant ontrols
are no longer appropriate. Preliminary omputer experiments with the lassial
Runge-Kuttamethod showthat orderof onvergene greaterthan two is attain-
able. But for these Runge-Kutta methods suitable seletion strategies have to
be studied in more detail. In this ontext, additional diÆulties arise if state
onstraints are present, beause these onstraints should be fullled also at the
intermediatestages ofthe Runge-Kuttasheme (as in[8℄).
FurtherresearhanbeondutedtowardsthestudyofRunge-Kuttashemes
asin[29℄, [13℄, [27℄,wherethe seletionstrategyismotivatedby multipleontrol
integrals. Inthe speialase oftwoseletionsperRunge-Kuttastep thisleads to
alternativeseletionsets oftype
u (i)
app (t
i +
1
h;^u);u (i)
app (t
i +
2 h;^u)
2
^
U U
U, whereUU orrespondstoase(iii)ofindependentseletionsinSetion3.2.
Thisset
^
U anbedesribedbynitelymanynonlinearinequalitiesandequalities,
whih an be easilyimposed as additionalonstraintsin the disretized optimal
ontrol problems.
Theproposedmethoditselfanbeeasilyadaptedtothe alulationof onvex
reahable sets for nonlinear dierentialinlusions. Forthe numerialsolution of
disretized optimal ontrol problems eÆient algorithms are available, f., e.g.,
[6℄, [16, 17℄. Inthe more generalase of nononvex reahable sets suitablemodi-
ationsofour approahhave tobestudied. Theoretialresultsinthis diretion
an be found in [12℄, [41℄, [40℄ for Runge-Kutta methods of order one and two.
A survey of other methods isgiven in[11℄ and [8℄.
However, those Runge-Kutta methods with appropriate seletion strategies,
whih show higher order of onvergene in the linear ase, are worth being in-
vestigated also in the nonlinear ase. In addition, these methods have to be
omparedwithset-valuedRunge-Kuttamethodsbasedonset arithmetis,f.[8℄,
whih work also on the general nonlinear ase. First steps in this diretion an
be found in[8, Example 5.3.1℄.
Referenes
[1℄ J.-P. Aubin and A. Cellina. Dierential Inlusions, volume 264 of
Grundlehren der mathematishen Wissenshaften. Springer Verlag, Berlin{
Heidelberg{New York{Tokyo, 1984.
[2℄ R. J. Aumann. Integrals of Set-Valued Funtions. J. Math. Anal. Appl.,
12(1):1{12, 1965.
ihbarer Mengen. Bayreuth. Math. Shr.,50:xxii + 248 S., 1995.
[4℄ R.Baierand F. Lempio. ComputingAumann's integral. In A.B. Kurzhan-
ski and V. M. Veliov, editors, Modeling Tehniques for Unertain Systems,
Proeedingsof a Conferenes heldin Sopron, Hungary,July6-10, 1992,vol-
ume18ofProgressinSystemsandControlTheory,pages71{92,Basel,1994.
Birkhauser.
[5℄ V. I. Blagodatskikh and A. F. Filippov. Dierentialinlusionsand optimal
ontrol. In Pro. Steklov Inst. Math., 4, pages 199{259. North-Holland,
Amsterdam, 1986.
[6℄ C. Buskens. Optimierungsmethoden und Sensitivitatsanalyse f ur optimale
Steuerprozesse mitSteuer- und Zustandsbeshrankungen. PhD thesis,Fah-
bereihMathematik,Westfalishe Wilhems-UniversitatMunster,1998.
[7℄ J.C. Buther. TheNumerialAnalysisof OrdinaryDierentialEquations|
Runge-Kutta and General Linear Methods. John Wiley and Sons,
Chihester{New York{Brisbane{Toronto{Singapore, 1987.
[8℄ I. A. Chahma. Set-valued disrete approximation of state-onstrained dif-
ferentialinlusions. Bayreuth. Math. Shr., 67:3{162, 2003.
[9℄ P. Diamond, P. Kloeden, A. Rubinov, and A. Vladimirov. Comparative
Properties of Three Metris in the Spae of Compat Convex Sets. Set-
Valued Anal., 5(3):267{289,1997.
[10℄ A.L.Donthev,W.W.Hager,andV.M.Veliov.Seond-OrderRunge-Kutta
ApproximationsinControlConstrainedOptimalControl. SIAM Journalon
NumerialAnalysis,38(1):202{226,2000.
[11℄ A.L.DonthevandF.Lempio.Dierenemethodsfordierentialinlusions:
A survey. SIAM Rev., 34(2):263{294,1992.
[12℄ A.L.DonthevandE.M.Farkhi.ErrorEstimatesforDisretizedDierential
Inlusions. Computing, 41:349{358, 1989.
[13℄ R. Ferretti. High-Order Approximations of Linear Control Systems via
Runge-KuttaShemes. Computing, 58(4):351{364,1997.
[14℄ A. F. Filippov. Dierential Equations with Disontinuous Righthand Sides.
Mathematis and Its Appliations (Soviet Series). Kluwer Aademi Pub-
lishers,Dordreht{Boston{London, 1988.
[15℄ J. E. Gayek. Approximating reahable sets for a lass of linear ontrolsys-
tems. Internat. J. Control, 43(2):441{453, 1986.
algebraishen Gleihungssystemen hoheren Indexes und ihre Anwendungen
in der Kraftfahrzeugsimulation und Mehanik. volume 61 of Bayreuther
Mathematishe Shriften, Bayreuth,2001.
[17℄ M. Gerdts. Diret Shooting Method for the Numerial Solution of Higher
IndexDAEOptimalControlProblems.Journal ofOptimizationTheoryand
Appliations, 117(2):267{294,2003.
[18℄ P. E. Gill,W.Murray, M. A.Saunders, and M.H. Wright. User's guide for
NPSOL 5.0: A FORTRAN pakage for nonlinear programming. Tehnial
Report NA 98-2,Department of Mathematis, Universityof California, San
Diego,California,1998.
[19℄ J.-B. Hiriart-Urruty and C. Lemarehal. Convex Analysis and Minimiza-
tion Algorithms I. Fundamentals, volume305 of Grundlehren der mathema-
tishen Wissenshaften. Springer, Berlin{Heidelberg{New York{London{
Paris{Tokyo{Hong Kong{Barelona{Budapest, 1993.
[20℄ J.E.Kelley. Theutting-planemethodforsolvingonvexprograms. J.So.
Ind. Appl. Math., 8:703{712,1960.
[21℄ K.C.Kiwiel.MethodsofDesentforNondierentiableOptimization,volume
1133ofLetureNotesinMath.Springer,Berlin-Heidelberg-NewYork.Tokyo,
1985.
[22℄ K.C. Kiwiel.A ConstraintLinearizationMethodforNondierentiableCon-
vex Minimization. Numer. Math.,51:395{414, 1987.
[23℄ A.B. KurzhanskiandI.Valyi. EllipsoidalCalulusforEstimationand Con-
trol. Systems & Control: Foundations& Appliations. Birkhauser, Boston{
Basel{Berline,1997.
[24℄ A. B. Kurzhanski and P. Varaiya. Ellipsoidal tehniques for reahability
analysis: internalapproximation. Systems Control Lett., 41:201{211, 2000.
[25℄ A. B. Kurzhanski and P. Varaiya. On Ellipsoidal Tehniques for Reaha-
bility Analysis. part I: External Approximations. Optim. Methods Softw.,
17(2):177{206,2002.
[26℄ A.B.KurzhanskiandP.Varaiya.OnEllipsoidalTehniquesforReahability
Analysis. part II: Internal Approximations Box-valued Constraints. Optim.
Methods Softw., 17(2):207{237,2002.
[27℄ P. E. Kloeden L. Grune. Higher order numerial shemes for aÆnely on-
trollednonlinear systems. Numer. Math., 89:669{690,2001.
Nonsmooth Optimization. In O. L. Mangasarian, R. R. Meyer, and S. M.
Robinson, editors, Nonlinear Programming 4, pages 245{282, New York,
1981.Aademi Press.
[29℄ F.LempioandV.Veliov.DisreteApproximationsofDierentialInlusions.
Bayreuth. Math. Shr., 54:149{232,1998.
[30℄ K.Malanowski,C. Buskens, and H.Maurer. Neessary ConditionsforOpti-
malControlProblemsInvolvingNonlinearDierentialAlgebraiEquations.
In Anthony Fiao, editor, Mathematial programming with data perturba-
tions, volume 195, pages 253{284. Dekker. Let. Notes Pure Appl. Math.,
1997.
[31℄ R. Miin. A modiation and an extension of Lemarehal's algorithm for
nonsmooth minimization. Math. Program. Study, 17:77{90, 1982.
[32℄ R. T. Rokafellar. Convex Analysis, volume 28 of Prineton Mathematial
Series.PrinetonUniversityPress,Prineton,NewYersey,2 nd
edition,1972.
[33℄ P. Saint-Pierre. Approximation of the viability kernel. Appl. Math. Optim.,
29:187{209,1994.
[34℄ K. Shittkowski. On the Convergene of a Sequential Quadrati Program-
ming Method with anAugmented Lagrangian Line Searh Funtion. Math.
Operationsforsh. Stat., Ser. Optimization,14(2):197{216,1983.
[35℄ K.Shittkowski. NLPQL:AFortransubroutineforsolvingonstrainednon-
linearprogramming problems. Ann. Oper. Res., 5:484{500,1985.
[36℄ H. Shramm. Eine Kombination von Bundle- und Trust-Region-Verfahren
zur Losung nihtdierenzierbarer Optimierungsprobleme. Bayreuth. Math.
Shr.,30,1989.
[37℄ L. M. Sonneborn and F. S. van Vlek. The bang-bang priniple for linear
ontrolproblems. SIAM J. Control, Ser. A, 2(2):151{159,1965.
[38℄ P. Varaiya. Reah set omputation using optimal ontrol. In M. K. Inan
and R. P. Kurshan, editors, Veriation of digital and hybrid systems. Pro-
eedings of the NATO ASI, Antalya, Turkey, May 26{June 6, 1997, volume
170of NATO ASI Ser,Ser. F,Comput. Syst.Si.,pages323{331.Springer,
2000.
[39℄ V.M.Veliov.Approximationstodierentialinlusionsbydisreteinlusions.
Working PaperWP-89-017, IIASA, Laxenburg, Austria, 1989.
ferentialinlusions. Systems Control Lett., 13:263{269,1989.
[41℄ V. M. Veliov. Seond Order Disrete Approximation to Linear Dierential
Inlusions. SIAM J. Numer. Anal., 29(2):439{451,1992.
[42℄ S.E.Wright.Primal-DualInterior-PointMethods.SIAM,Philadelphia,PA,
1997.
[43℄ G.Zoutendijk.MethodsofFeasibleDiretions.ElsevierPublishingCompany,
Amsterdam, Netherlands, 1960.
