Vorabdruck!

Numerische Losung eines optimalen Steuerproblemsaus der chemischen Reaktionskinetik

mit der NAG-Bibliothek, MATLAB und Mathematica

Siegfried Dietze, TU Dresden, Institut fur Numerische Mathematik, Marz 2005

Inhaltsubersicht: Im vorliegenden Skript wird ein einfaches optimales Steuerproblem ausder chemischen Reaktionskinetik mit der NAG-Fortran-Bibliothek, Matlab und Mathema-tica numerisch gelost.Das Steuerproblem hat einen endlichdimensionalen Steuerraum, und Steuerbeschrankun-gen sowie andere Beschrankungen werden nicht einbezogen. Demnach handelt es sich umeine endlichdimensionale Optimierungsaufgabe ohne Nebenbedingungen, so daß das Steu-erproblem mit Software zur freien Minimierung gelost werden kann. Eine Besonderheit ist,daß zur Berechnung jedes Zielfunktionswertes eine Anfangswertaufgabe bei gewohnlichenDifferentialgleichungen gelost werden muß.Die Losung der Optimierungsaufgabe mit den drei verwendeten Softwaresystemen wirdausfuhrlich erortert, so daß das Skript - eventuell noch erganzt durch das Kapitel Softwareaus dem Lehr- und Ubungsbuch zur Numerischen Mathematik [8] - eine gute Vorlage zurLosung eigener anspruchsvollerer Aufgaben mit einem der drei Softwaresysteme darstellt.Das erwahnte Kapitel aus [8] stellt eine vergleichende Ubersicht uber Standardsoftware zurNumerik dar und enthalt kurze aber aussagekraftige Einfuhrungen zur NAG-Bibliothek,zu Matlab und Mathematica.Durch die Bearbeitung ein- und derselben Aufgabe mit drei Softwaresystemen werden dieim Grunde bekannten Vor- und Nachteile der Software besonders ersichtlich. Das erleich-tert die passende Softwareauswahl bei der Losung eigener Aufgaben.

Inhalt:

1. Die Optimierungsaufgabe: Maximale Ausbeute bei einer chemischen Reaktion

2. Bearbeitung der Optimierungsaufgabe mit der NAG-Bibliothek

3. Bearbeitung der Optimierungsaufgabe mit Matlab

4. Bearbeitung der Optimierungsaufgabe mit Mathematica

5. Rechenzeiten

6. Schlußfolgerungen

Literatur

1

1 Die Optimierungsaufgabe: Maximale Ausbeute bei

einer chemischen Reaktion

Bei der Optimierungsaufgabe handelt es sich um das einfachste Beispiel aus der chemischenReaktionskinetik. Betrachtet wird die Folgereaktion

A1k1−→ A2

k2−→ A3

der Stoffe A1, A2, A3 mit den von der Temperatur abhangigen Geschwindigkeitskonstantenk1, k2 (vgl. [4], [5], [6]). Die Reaktion erfolgt in einem Rohrreaktor der Lange L. Die StoffeA1, A2, A3 sind gasformig und durchstromen den Reaktor von links nach rechts mit derkonstanten Geschwindigkeit v, gleichmaßig uber den Rohrquerschnitt verteilt.

- v

- l0 L

- t0 t1 = L/v

l

t

Als unabhangige Variable wird die Zeit t = l/v verwendet. Die Temperatur im Reaktorzur Zeit t, d.h. an der t entsprechenden Stelle l, sei T(t), gemessen in Kelvin und uber denRohrquerschnitt konstant. Fur die Geschwindigkeitskonstanten ki gelte

ki = ki(t) = ki0e−Ei

RT (t) (i = 1, 2). (1)

Hier ist R die Gaskonstante und Ei eine positive Konstante. Strebt die Temperatur gegenNull, so strebt die Geschwindigkeitskonstante ki gegen Null. Wachst die Temperatur, sowachst ki und strebt mit t→ +∞ gegen ki0. Es sei nun

xi = xi(t) = ��Ai �� (i = 1, 2, 3)

die Konzentration des Stoffes Ai zur Zeit t. Auf der Basis des Massenwirkungsgesetzes kannman zeigen, daß die Konzentrationen xi(t) dem Anfangswertproblem

x1 = −k1(t)x1 , x1(0) = x10

x2 = k1(t)x1 − k2(t)x2, x2(0) = x20 (0 ≤ t ≤ t1) (2)

x3 = k2(t)x2, x3(0) = x30

genugen. Der Literatur folgend transformieren wir noch das Zeitintervall [0, t1] auf dasIntervall [0, 1], was jedoch nicht zwingend ist: Wenn wir dabei fur die neu eingefuhrtenGroßen t := t/t1, T (t) := T (t1t) = T (t) und xi(t) := xi(t1t) = xi(t) wieder die alten

2

Bezeichnungen benutzen, also die Schlange ˜ uberall weglassen, gehen die beiden erstenGleichungen des Anfangswertproblems (2) uber in

x1 = −t1k1(t)x1 , x1(0) = x10

x2 = t1k1(t)x1 − t1k2(t)x2, x2(0) = x20(0 ≤ t ≤ 1). (3)

Die dritte Gleichung von (2) spielt im weiteren keine Rolle, weil die Ausbeute von A2

interessiert, d.h. x2(1).Die Ausbeute x2(1) = x2(1;T (·)) von A2 hangt ab von der gewahlten TemperaturfunktionT (t), 0 ≤ t ≤ 1. Dies fuhrt auf die Optimierungsaufgabe

−x2(1;T (·)) −→ minT (·)∈D

. (4)

Dabei ist D die Menge der folgenden auf [0, 1] erklarten Treppenfunktionen: Es sei neine fest gewahlte naturliche Zahl. Dadurch sind die Schrittweite und die aquidistantenGitterpunkte

h = 1/n, tj = jh (j = 0, . . . , n)

festgelegt. Eine Treppenfunktion T (t) = T (t;u), 0 ≤ t ≤ 1 ist dann durch einen Vektoru = (u1, . . . , un)T ∈ IRn durch

T (t) = T (t;u) =

{

uj fur t ∈ [tj−1, tj) (j = 1, . . . , n− 1)un fur t ∈ [tn−1, tn]

gegeben. Das folgende Bild zeigt eine Treppenfunktion im Fall n = 4.

- t

6T

t0 = 0 t1 t2 t3 t4 = 1

u1

u2

u3

u4

Wir weisen noch einmal darauf hin, daß die naturliche Zahl n fest gewahlt wird. Folglichsind auch die Gitterpunkte tj fest. Die Menge D der Temperatur-Treppenfunktionen ergibtsich durch Variation des Vektors u und die Optimierungsaufgabe (4) kann in der Form

−x2(1;T (·;u)) −→ minu∈IRn

. (5)

3

geschrieben werden.Bei (5),(3),(1) handelt es sich um eine endlichdimensionale Optimierungsaufgabe ohne Ne-benbedingungen. Zur Berechnung eines Zielfunktionswertes (5) zum Vektor u ∈ IRn ist dieAnfangswertaufgabe (3) zu losen, wobei die Koeffizienten ki(t) (1) entsprechend der Trep-penfunktion T (·) = T (·;u) zu wahlen sind. Wir werden Minimierungs-Software verwenden,bei der der Nutzer keine Ableitungen bereitstellen muß und bei der erforderliche Ableitun-gen intern uber Differenzenquotienten berechnet werden. Die ersten partiellen Ableitungender Zielfunktion konnten uber das zu (3) adjungierte Differentialgleichungssystem berech-net werden (vgl. [5], [11]), was bei großerem n zu Effektivitatssteigerungen fuhrt. Daraufgehen wir hier nicht ein.Die dem Vektor u ∈ IRn entsprechende Anfangswertaufgabe (3) werden wir numerischlosen. Dabei ist zu beachten, daß die rechte Seite des Differentialgleichungssystems in tstuckweise konstant ist mit den Sprungstellen tj. Die Losung der Anfangswertaufgabe istzwar auf ganz [t0, tn] = [0, 1] stetig, jedoch an den Stellen tj nicht differenzierbar. Manlost die Anfangswertaufgabe, indem man sie zunachst auf dem Intervall [t0, t1] lost, wichtigsind die Endwerte x1(t1), x2(t1). Anschließend lost man die Differentialgleichung (3) aufdem Intervall [t1, t2], wobei die Endwerte vom ersten Intervall als Anfangswerte ubernom-men werden, usw. Dieses Anstuckeln der Losung ist auch beim Einsatz der numerischenVerfahren zu realisieren. Denn sonst ist nicht sicher, daß das numerische Verfahren mitder zugehorigen Ordnung konvergiert. Die dem Vektor u ∈ IRn entsprechende Anfangs-wertaufgabe (3) ist stuckweise linear mit konstanten Koeffizienten und konnte daher auchgeschlossen gelost werden - wiederum durch Anstuckeln der Losung, sofern symbolischesRechnen in der verwendeten Software moglich ist. Von dieser Moglichkeit haben wir vorallem deshalb abgesehen, weil das Differentialgleichungssystem fur die Konzentrationen deran einer chemischen Reaktion beteiligten Stoffe im allgemeinen nichtlinear ist, aber auchaus Rechenzeitgrunden.

Bemerkung 1: Die Aufgabe (5),(3),(1) bzw. (4),(3),(1) ist ein einfaches optimales Steu-erproblem ohne Steuer- und sonstige Beschrankungen. Dieses kann als Diskretisierung desSteuerproblems aufgefaßt werden, das aus (4),(3),(1) entsteht, wenn dort an Stelle der Trep-penfunktionen allgemeinere Funktionen zugelassen werden. Wir erwahnen, daß zusatzlichauch die Zustandsfunktion x(·) = (x1(·), x2(·))T durch Anwendung eines Diskretisierungs-verfahrens auf die Anfangswertaufgabe (3) diskretisiert werden konnte (vgl. [1], [3], [11]).Wendet man auf die wie auch immer diskretisierte Aufgabe einen Loser fur endlichdimensio-nale Optimierungsaufgaben an, so spricht man von einem direkten Verfahren zur Losung desSteuerproblems. Neben den direkten Verfahren gibt es indirekte Verfahren. Ausgangspunktdafur ist die bekannte Optimalitatsbedingung fur Steuerprobleme, das Maximumprinzip[7]. Beim vorliegenden Steuerproblem in nichtdiskretisierter Form fuhrt das Maximumprin-zip auf eine nichtlineare Randwertaufgabe bei gewohnlichen Differentialgleichungen (vgl.[4], [6]), die z.B. mit dem Schießverfahren [10] numerisch gelost werden kann.

Bemerkung 2: Dem Verfasser ist nicht bekannt, ob die Optimierungsaufgabe (5),(3),(1)losbar ist. Dasselbe gilt fur das zugehorige nichtdiskretisierte Steuerproblem aus Bemer-kung 1. Losbarkeitsaussagen siehe z.B. [2]. Vom praktischen Problem her ist allerdings klar,

4

daß die Zielfunktion in (5) wenigstens nach unten beschrankt ist.

Fur unsere Rechnungen verwenden wir folgende Daten aus der Literatur: x10 = 0.53 mol dm−3,x20 = 0.43 mol dm−3, k10 = 0.892 · 109 s−1, k20 = 0.768 · 1016 s−1, E1 = 18000 cal mol−1,E2 = 30000 cal mol−1, R = 2 cal mol−1 K−1, t1 = 480 s.

2 Bearbeitung der Optimierungsaufgabe mit der NAG-

Bibliothek

Wir verwenden die fl90, d.h. die Fortran-90-Bibliothek der Firma NAG, in der Version mitdoppelter Genauigkeit. Eine allgemeine Einfuhrung in diese Bibliothek soll hier nicht ge-ben werden. Dazu verweisen wir auf die ausfuhrliche Dokumentation, die auch im Internetunter http://www.nag.co.uk eingesehen werden kann. Eine kurze Einfuhrung einschließlichvollstandiger Bearbeitung eines Beispiels zur Nutzung der Bibliothek auf einer Worksta-tion IBM RS/6000 mit dem Compiler AIX XL von IBM (! Endung .f und nicht .f90 furden Quelltext) findet man in [8]. Auf einer solchen Maschine wurden auch die in die-sem Abschnitt angegebenen Ergebnisse erzielt. Wie im Abschnitt 1 diskutiert besteht diezu losende Optimierungsaufgabe im wesentlichen aus den Teilaufgaben Minimierung undAnfangswertaufgaben bei gewohnlichen Differentialgleichungen . Beide Teilaufgaben sindin der fl90 vertreten. Zur Losung der Optimierungsaufgabe ist ein Fortran-90-Programmerforderlich, in dem die fur die genannten Teilaufgaben entsprechenden Unterprogrammeaus der fl90 gerufen werden. Bei der Erstellung dieses Programmes kann auf die von NAGmitgelieferten Beispiel-Hauptprogramme fur die einzelnen Bibliotheksbausteine zuruckge-griffen werden, so daß nur geringe Kenntnisse in Fortran erforderlich sind. Im weiterenwerden das Programm zur Losung der Optimierungsaufgabe diskutiert und die damit er-zielten Ergebnisse angegeben. Auf die Abarbeitung des Programmes gehen wir nur amRande ein.Den Teilaufgaben Minimierung und Anfangswertaufgaben bei gewohnlichen Differential-gleichungen entsprechend interessieren von der fl90 der Teil 9: Optimierung und der Teil12: Gewohnliche Differentialgleichungen.

a) Wahl des Unterprogrammes zur Minimierung: Der Teil 9: Optimierung bestehtzur Zeit (Release 4) aus sechs Modulen. Fur uns ist Module 9.3: nag nlp - Nonlinear Pro-gramming zur Losung nichtlinearer Optimierungsaufgaben ohne und mit Nebenbedingun-gen wichtig. Dieser Modul enthalt die eigentliche Optimierungsprozedur nag nlp sol, denabgeleiteten Datentyp nag nlp cntrl wp fur Steuerparameter von nag nlp sol (z.B. Parame-ter fur Ausgabe, Genauigkeitsschranken) sowie die Initialisierungsprozedur nag nlp cntrl initfur nag nlp cntrl wp. Uns interessiert nur nag nlp sol, weil wir mit voreingestellten Steuer-parametern arbeiten werden. Bei nag nlp sol handelt es sich um die Implementierung einesSQP-Verfahren, siehe Dokumentation und z.B. [1]. SQP-Verfahren sind Iterationsverfahrenvom Typ

neue Iterierte = alte Iterierte + Schrittweite · Suchrichtung,wobei in jedem Schritt zur Bestimmung der Suchrichtung ein quadratisches Optimierungs-

5

problem zu losen ist. Das Beispiel-Hauptprogramm nag nlp ex05 von NAG fur die Mi-nimierung einer konkreten Funktion ohne Nebenbedingungen mittels nag nlp sol ist einegute Grundlage fur das von uns zu erstellende Programm. Dieses Beispiel-Hauptprogrammist in der Dokumentation nicht vollstandig enthalten , es ist aber als Datei nag nlp ex05.fBestandteil der NAG-Bibliothek und sollte eingesehen werden. Wir verzichten hier auf dieWiedergabe dieses Programmes, wollen aber auf den Aufruf von nag nlp sol in der fur unsverwendeten Form und die Bedeutung der Argumente etwas naher eingehen, ausfuhrlicheErlauterungen siehe Dokumentation:

CALL nag nlp sol(obj fun,u,obj f,obj deriv=.false.)

obj fun Unterprogramm fur Zielfunktionswert und Gradient:

subroutine obj fun(first call,u,finish,obj f,obj grad)first call fur uns unwichtigu Stelle, wo Funktionswert und Gradient berechnet werden sollenfinish fur uns unwichtigobj f Funktionswert an der Stelle uobj grad Gradient an der Stelle u

Wegen obj deriv=.false. im rufenden Programm wird der Gradient automa-tisch mittels Differenzenquotienten approximiert.(Bei obj deriv=.true. – dies entspricht der Voreinstellung – mußte dieAbleitung in obj fun bereitgestellt werden, vgl. Beispielhauptprogrammnag nlp ex05.)

u Startvektor bzw. Losungobj f Funktionswert zur Losungobj deriv optionales Argument; spezifiziert, ob Gradient vom Nutzer in obj fun be-

reitgestellt wird oder dort automatisch mittels Differenzenquotienten ap-proximiert werden soll.

Es gibt zahlreiche weitere optionale Argumente zu Nebenbedingungen des Optimierungs-problems und zu Verfahrensparametern. Diese spielen fur uns keine Rolle, weil wir mit Vor-einstellungen arbeiten. Insbesondere ist keine Nebenbedingungen voreingestellt. Die Vorein-stellungen der optionalen Argumente wie auch die Voreinstellungen der in nag nlp cntrl wpzusammengefaßten Steuerparameter findet man in der Dokumentation. Einige sind aberauch aus den unten ausgedruckten Ergebnissen indirekt ablesbar.

b) Wahl des Unterprogrammes zur Losung der Anfangswertaufgabe: Der Teil12 der fl90: Gewohnliche Differentialgleichungen besteht zur Zeit nur aus Module 12.1:nag ivp ode rk - Solution of InitialValue Problems for ODE’s by Runge-Kutta Methods.Dazu ist zu bemerken, daß die fl90 noch im Aufbau begriffen ist. In der umfangreichenFortran-77-Bibliothek von NAG ist wesentlich mehr uber gewohnliche Differentialgleichun-gen enthalten, auch Loser fur steife Differentialgleichungen und Randwertaufgaben. BeiBedarf kann eine Routine aus dieser Bibliothek in ein Fortran-90-Programm eingebundenwerden, vgl. Dokumentation und ein Beispiel auf Workstation IBM RS/6000 mit Compiler

6

AIX XL von IBM in [8]. Da unsere Anfangswertaufgabe nicht steif ist, reicht der Modul12.1. Dieser enthalt die beiden Prozeduren nag rk setup und nag rk interval zur Losungder Anfangswertaufgabe, den abgeleiteten Datentyp nag rk comm wp zur Kommunikati-on zwischen den Prozeduren sowie weitere hier uninteressante Prozeduren. Die Prozedurnag rk setup ist vor dem eigentlichen Anfangswertloser nag rk interval zu rufen und dientder Initialisierung (vergleiche unten angegebene Aufrufe und Parameterlisten). Wie bereitsder Moduluberschrift zu entnehmen, handelt es sich bei nag rk interval um eine Implemen-tierung von Runge-Kutta-Verfahren. Die Schrittweite wird nach dem Einbettungsprinzipgesteuert, siehe z.B. [9], [10]. Es sind die Ordnungen 2(3), 4(5) und 7(8) wahlbar, wir arbei-ten mit der Voreinstellung 4(5). Als Grundlage fur den Einsatz des Anfangswertlosers furunsere Belange verweisen wir auf das Beispiel-Hauptprogramm nag ivp ode rk ex01 vonNAG, das eingesehen werden sollte. Auf die Wiedergabe dieses Programmes verzichtenwir. Es folgen die Aufrufe von nag rk setup und nag rk interval fur die Anfangswertaufga-be y = f(t, y), y(t0) = y0 in der von uns verwendeten Form und eine Kurzbeschreibungder Argumente, ausfuhrliche Erlauterungen siehe Dokumentation:

CALL nag rk setup(t start,y start,t end,tol,thresh,comm)fur Bereitstellung von Problemdaten, Daten fur Fehlerkontrolle, . . .

t start Anfangszeit t0y start Anfangswert y0

t end Endzeit tetol fur Fehlerkontrollethresh Vektor von thresholds=Schwellenwerten fur Fehlerkontrolle

Es wird die Schrittweite uber den lokalen Fehler nach dem Einbettungs-prinzip gesteuert. Fur die Komponente i des lokalen Fehlers wird die Feh-lerschranke tol∗max{ηi, thresh(i)} benutzt. Dabei ist ηi ein Naherungswertfur die Losungskomponente yi(t) an der jeweiligen Stelle. Es wird also relativgetestet fur ηi ≥ tresh(i), sonst absolut.

comm ist vom Typ nag rk comm wp und dient der Kommunikation zwischen denProzeduren, wird intern belegt.

CALL nag rk interval(f,t want,t got,y got,yp got,comm)

f Funktion, rechte Seite der Differentialgleichungt want gewunschtes t zwischen t0 und Endzeit te, wo Losungswert gesuchtt got (Ausgabewert,) erhaltenes t. Im Normalfall gilt t got=t want. Es kann aber

sein, daß t want nicht erreicht wurde.y got Losung an der Stelle t gotyp got Ableitung der Losung an der Stelle t gotcomm wie oben

Beide Prozeduren haben auch optionale Argumente, z.B. hat nag rk setup ein optionalesArgument zur Wahl der Methode.

c) Programm-Quelltext reaktor.f: Wir geben jetzt den Programm-Quelltext zur Losungder Optimierungsaufgabe aus Abschnitt 1 an. Wie bereits oben erwahnt, ist das Programm

7

auf der Grundlage der Beispielhauptprogramme nag nlp ex05 und nag ivp ode rk ex01 vonNAG erstellt worden, wobei die dort enthaltenen Kommentare nicht konsequent fortgefuhrtworden sind. Der Quelltext ist im File reaktor.f gespeichert:

MODULE reaktor_mod

! .. Implicit None Statement ..

IMPLICIT NONE

! .. Intrinsic Functions ..

INTRINSIC KIND

! .. Parameters ..

INTEGER, PARAMETER :: wp = KIND(1.0D0)

INTEGER :: n

REAL (wp) :: k10,k20,k1,k2

REAL (wp) :: t(51)

COMMON n, /ber1/ t,k10,k20

CONTAINS

SUBROUTINE obj_fun(first_call,u,finish,obj_f,obj_grad)

! .. Use Statements ..

USE nag_ivp_ode_rk, ONLY : nag_rk_comm_wp => nag_rk_comm_dp, &

nag_rk_setup, nag_rk_interval, nag_deallocate

! .. Implicit None Statement ..

IMPLICIT NONE

! .. Scalar Arguments ..

REAL (wp), INTENT (OUT) :: obj_f

LOGICAL, INTENT (INOUT) :: finish

LOGICAL, INTENT (IN) :: first_call

! .. Array Arguments ..

REAL (wp), OPTIONAL, INTENT (INOUT) :: obj_grad(:)

REAL (wp), INTENT (IN) :: u(:)

INTEGER, PARAMETER :: m = 2

INTEGER :: j

REAL (wp) :: t_start, t_end, tol, t_got

TYPE (nag_rk_comm_wp) :: comm

REAL (wp) :: thresh(m), yp_got(m), y_got(m), y_start(m)

! .. Executable Statements ..

tol = 1.0E-7_wp

thresh = 1.0E-10_wp

y_got = (/ 0.53_wp, 0.43_wp /)

8

DO j = 1, n

k1 = k10*EXP(-9000.0_wp/u(j))

k2 = k20*EXP(-15000.0_wp/u(j))

t_start = t(j)

t_end = t(j+1)

y_start = y_got

CALL nag_rk_setup(t_start,y_start,t_end,tol,thresh,comm)

CALL nag_rk_interval(f,t_end,t_got,y_got,yp_got,comm)

END DO

obj_f = -y_got(2)

CALL nag_deallocate(comm) ! Free structure allocated by NAG fl90

END SUBROUTINE obj_fun

FUNCTION f(t,y)

! .. Implicit None Statement ..

IMPLICIT NONE

! .. Intrinsic Functions ..

INTRINSIC SIZE

! .. Scalar Arguments ..

REAL (wp), INTENT (IN) :: t

! .. Array Arguments ..

REAL (wp), INTENT (IN) :: y(:)

! .. Function Return Value ..

REAL (wp) :: f(SIZE(y))

! .. Executable Statements ..

f(1) = -k1*y(1)

f(2) = k1*y(1)-k2*y(2)

END FUNCTION f

END MODULE reaktor_mod

PROGRAM reaktor

! .. Use Statements ..

USE nag_examples_io, ONLY : nag_std_out

USE nag_nlp, ONLY : nag_nlp_sol

USE reaktor_mod, ONLY : obj_fun, wp

! .. Implicit None Statement ..

IMPLICIT NONE

INTEGER j,n

REAL (wp) :: h, obj_f, k10, k20

REAL (wp), ALLOCATABLE :: u(:)

REAL (wp) :: t(51)

9

COMMON n, /ber1/ t,k10, k20

! .. Executable Statements ..

WRITE (*,*) ’n=’

READ (*,*) n

WRITE (*,*) n

ALLOCATE(u(n))

h=1.0_wp/n

DO j = 1, n+1

t(j) = (j-1)*h

END DO

k10 = 480.0_wp*0.892E9_wp

k20 = 480.0_wp*0.768E16_wp

WRITE (nag_std_out,*) ’Example Program Results for reaktor’

u=350.0_wp

! Solve the problem

CALL nag_nlp_sol(obj_fun,u,obj_f,obj_deriv=.false.)

!Erzeugung einer Datei u.d mit Inhalt n und u

OPEN(10,FILE=’u.d’)

WRITE(10,*) n

WRITE(10,*) u

END PROGRAM reaktor

Das Programm besteht aus dem Modul reaktor mod und dem Hauptprogramm reaktor. ImModul werden Funktionswert und Gradient der zu minimierenden Funktion (SUBROUTI-NE obj fun), die rechte Seite des Differentialgleichungssystems (FUNCTION f) und derTypparameterwert wp fur doppelte Genauigkeit bereitgestellt, und im Hauptprogrammwird im wesentlichen der Minimierer nag nlp sol gerufen. In beiden Programmeinheitensind die gemeinsamen Speicherbereiche n (unbenannter COMMON-Bereich) und /ber1/t,k10,k20 (benannter COMMON-Bereich) vereinbart. Die entsprechenden Variablen wer-den im Hauptprogramm gesetzt und sind dann im Modul verfugbar. Die Trennung in diezwei Bereiche wurde vorgenommen, weil sonst beim Ubersetzen eine Effektivitatswarnungerfolgte, vermutlich wegen der Mischung von INTEGER- und REAL-Variablen.Im Modul sind Vereinbarungen global gultig, also auch in den Untereinheiten. Die beidenAufrufe fur den Anfangswertloser CALL nag rk setup(t start,y start,t end,tol,thresh,comm)und CALL nag rk interval(f,t end,t got,y got,yp got,comm) haben wir bereits in b) ausfuhr-lich diskutiert. Mit der j-Schleife um den Anfangswertloser wird das im Abschnitt 1 be-schriebene Anstuckeln der Losung realisiert. Dabei wird stillschweigend angenommen, daßes keine Probleme bei der Losung auf dem Intervall [t start,t end]=[t(j), t(j+1)] gibt, d.h.,es wird t got=t end angenommen. Dies mußte eigentlich uberpruft werden.Im Hauptprogramm besteht eine gewisse Diskrepanz bei der Vereinbarung von u als dy-namisches Feld und t als statisches Feld. Weil t zu einem COMMON-Bereich gehort, kannes nicht dynamisch vereinbart werden. Aus der Vereinbarung von t ist ubrigens fur die

10

Stufenzahl n der Treppenfunktionen die Bedingung n≤50 ersichtlich, was im Programmnicht uberpruft wird. Die Zahl n wird eingelesen. Dagegen ist der Startvektor u fur den Mi-nimierer fest eingetragen. Der fur alle u-Komponenten gleiche Wert 350 hat sich bewahrt.Den Aufruf fur den Minimierer CALL nag nlp sol(obj fun,u,obj f,obj deriv=.false.) habenwir bereits in a) ausfuhrlich diskutiert. Neben der voreingestellten Ausgabe des Minimie-rers werden n und u, d.h. bei erfolgreicher Arbeit des Minimierers die ermittelte Nahe-rungslosung, in die Datei u.d geschrieben. Auf diese Datei konnte z.B. von einem Graphik-system aus bei der graphischen Darstellung der optimalen Temperatur-Treppenfunktionund/oder der zugehorigen Konzentrationen x1(t), x2(t) zuruckgegriffen werden. Darauf ge-hen wir nicht weiter ein. Bei der Bearbeitung der Optimierungsaufgabe mittels Matlab undMathematica werden wir die graphische Darstellung “der Losung” ausfuhrlich behandeln.

d) Ubersetzen und Ausfuhren von reaktor.f: Auf der Workstation IBM RS/6000steht der Fortran-90-Compiler AIX XL mit integriertem Linker von IBM unter dem Na-men xlf90 zur Verfugung. Das aktuelle Verzeichnis beinhaltet zunachst nur das File re-aktor.f, vgl. erste Eingabeaufforderung dietze@rcs8: des folgenden Bildschirmausschnittes.Das Unix-Kommando ls veranlaßt die Ausgabe des Inhalts des aktuellen Verzeichnisses.Der Compileraufruf fur reaktor.f steht nach der zweiten Eingabeaufforderung. Er beginntmit xlf90. Es folgen der nach -I angegebene Pfad zu den Modul-Symbol-Files der fl90 (dasabschließende backslash bewirkt die Fortsetzung der Eingabe auf der nachsten Zeile), dermit -L notierte Pfad zur Bibliothek libnagfl90.a, welche die ubersetzten Programme derfl90 enthalt, die nach -l angefuhrte Bibliothek, wobei vom Bibliotheksnamen die ersten dreiZeichen lib und die beiden letzten Zeichen .a weggelassen werden, und schließlich das zu be-arbeitende Programm reaktor.f. Nach Bestatigung des Kommandos erscheinen zeitversetztdie im Anschluß an das Kommando folgenden Compilermeldungen. Das aktuelle Verzeich-nis enthalt jetzt neben dem Programm reaktor.f das zugehorige Maschinenprogramm mitdem voreingestellten Namen a.out und das Modul-Symbol-File reaktor mod.mod. Fur dieAbarbeitung von a.out werden die ubrigen Files nicht mehr benotigt.

dietze@rcs8: ls

reaktor.f

dietze@rcs8: xlf90 -I/licsoft/nag-fl90/nag_mod_dir\

> -L/licsoft/nag-fl90 -lnagfl90 reaktor.f

** reaktor_mod === End of Compilation 1 ===

** reaktor === End of Compilation 2 ===

1501-510 Compilation successful for file reaktor.f.

dietze@rcs8: ls

a.out* reaktor.f reaktor_mod.mod

Der folgende Bildschirm zeigt den Aufruf von a.out und die erzielten Ergebnisse im Falln = 4. Nach der Eingabe von 4 fur n wird die Uberschrift Example Program Results forreaktor ausgegeben. Es folgt der voreingestellte recht umfangreiche Ausgabeteil des Mini-mierers mit einer Liste der Parameter, Informationen uber Differenzenquotienten der Ziel-funktion im Startpunkt, Aussagen uber die Iterierten des SQP-Verfahrens und schließlich

11

mit der ermittelten Naherungslosung. In jedem Iterationsschritt des SQP-Verfahrens wirdfur die Losung des quadratischen Teilproblems nur eine minor iteration benotigt, weil dieOptimierungsaufgabe keine Nebenbedingungen hat. Der Parameter optim tol= 5.36E-12ist entscheidend fur die erreichte Genauigkeit, vgl. Dokumentation. Es gilt u.a. fur die vor-letzte Iterierte u und die letzte Iterierte u, d.h. die Naherungslosung,‖ u− u ‖≤

√optim tol(1+ ‖ u ‖).

dietze@rcs8: a.out

n=

4

Example Program Results for reaktor

Parameters

----------

list................... .true. lt80_char.............. .true.

unit................... 6

linear constraints..... 0 variables.............. 4

nonlinear constraints.. 0

inf_bound.............. 1.00E+20 cold_start............. .true.

inf_step............... 1.00E+20 eps (machine precision) 2.22E-16

step_limit............. 2.00E+00 hessian................ .true.

major_print_level...... 10 minor_print_level...... 0

major_iter_lim......... 50 minor_iter_lim......... 50

linesearch_tol......... 9.00E-01 fun_prec............... 8.16E-15

optim_tol.............. 5.36E-12 crash_tol.............. 1.00E-02

obj_deriv.............. .false. con_deriv.............. .true.

cheap_test............. .true. initial_x.............. .false.

Difference intervals to be computed.

The user sets 0 out of 4 objective gradient elements.

Each iteration, 4 gradient elements will be estimated numerically.

Computation of the finite-difference intervals

----------------------------------------------

12

J X(J) Forward DX(J) Central DX(J) Error est.

1 3.50E+02 7.392725E-06 6.342131E-04 6.382255E-09

2 3.50E+02 6.930966E-06 6.342131E-04 6.807458E-09

3 3.50E+02 6.690444E-06 6.342131E-04 7.052187E-09

4 3.50E+02 6.556590E-06 6.342131E-04 7.196158E-09

Maj Mnr Step Objective Norm Gz Cond Hz

0 1 0.0E+00 -4.451819E-01 1.9E-02 1.0E+00

1 1 1.6E+03 -6.718635E-01 2.5E-03 1.0E+00

2 1 1.0E+00 -6.782556E-01 1.1E-03 1.1E+00

3 1 1.0E+00 -6.794056E-01 2.9E-04 1.1E+00

4 1 1.0E+00 -6.794826E-01 4.5E-05 1.1E+00

5 1 1.0E+00 -6.794850E-01 3.1E-05 2.0E+00

6 1 1.0E+00 -6.794875E-01 8.5E-06 2.0E+00

7 1 1.0E+00 -6.794876E-01 2.0E-06 1.8E+00

8 1 1.0E+00 -6.794876E-01 5.9E-08 1.8E+00

9 1 1.0E+00 -6.794876E-01 2.2E-08 2.4E+00

Exit from nag_nlp_sol after 9 major iterations,

10 minor iterations.

Varbl State Value Lower Bound Upper Bound Lagr Mult Slack

V 1 FR 339.487 None None

V 2 FR 337.811 None None

V 3 FR 336.822 None None

V 4 FR 336.162 None None

Exit nag_nlp_sol - Optimal solution found.

Final objective value = -0.6794876

Nach der Abarbeitung von a.out enthalt das aktuelle Verzeichnis zusatzlich zu den bereitsvorhandenen Files das im Programm erzeugte File u.d. Dieses hat im Fall n = 4 denfolgenden Inhalt:

4

339.486835758609971 337.811156536562123 336.822088838796390 336.161800689446864

13

3 Bearbeitung der Optimierungsaufgabe mit Matlab

Eine kurze Einfuhrung in Matlab findet man z.B. in [8]. Die Arbeit mit Matlab wird durchdie integrierte umfangreiche Hilfe und Dokumentation wesentlich unterstutzt.

a) Wahl des Minimierers: Als Minimierer verwenden wir das Matlab-Programmfminunc aus der Matlab-Optimization Toolbox, die nicht zum Basispaket von Matlabgehort. Der von uns verwendete Aufruf lautet

[uopt,fopt]=fminunc(’fwert’,u0,options)- vgl. script file reaktor.m unten in c). Das Ergebnis [uopt,fopt] enthalt bei erfolgreicher Mi-nimierung die Minimumstelle uopt und den zugehorigen Funktionswert fopt. Der Parameter’fwert’ steht fur das function file fwert.m, in dem der Funktionswert der zu minimierendenFunktion bereitgestellt wird, und u0 ist der Startvektor des Minimierungsverfahrens. DerParameter options enthalt in unserem Fall Genauigkeitsschranken und legt die Minimie-rungsmethode und die Art der Ausgabe fest. Options wird vor dem Aufruf von fminuncmit Hilfe des Matlab-Programmes optimset gesetzt, s. auch die entsprechenden Kommen-tare in reaktor.m. So wird z.B. durch ’LargeScale’ und ’off’ ein Quasi-Newton-Verfahren(vgl. [1]) gewahlt, das zur Losung von Optimierungsaufgaben ohne Nebenbedingungen biszu mittlerer Dimension geeignet ist und fur das der Nutzer keine Ableitungen bereitstellenmuß. Der Minimierer fminunc verwendet weitere Steuergroßen. Diese tauchen nicht explizitauf, weil wir mit Voreinstellungen arbeiten. Die vollstandige Dokumentation zu fminuncund optimset findet man in der Hilfe zu Matlab.Wir erwahnen noch, daß wir bewußt auf den Einsatz des Minimierers fminsearch aus demBasispaket von Matlab verzichtet haben. Dort wird das Suchverfahren von Nelder/Mead,das sogenannte Simplexverfahren, verwendet. Dieses Verfahren, das nur Funktionswertebenutzt, ist zwar bekanntermaßen recht robust aber fur Aufgaben großerer Dimension ausRechenzeitgrunden nicht geeignet [1]. Im Falle unserer Optimierungsaufgabe fuhrt der Ein-satz von fminsearch an Stelle von fminunc bei Verwendung gleicher Genauigkeitsschrankenim Fall der Dimension n = 5 zur etwa vierfachen und bei n = 10 bereits zur zehnfachenRechenzeit, siehe die in Abschnitt 5 angegebenen Rechenzeiten.

b) Wahl des Losers fur die Anfangswertaufgabe: Im Basispaket von Matlab sind ver-schiedene Loser fur nichtsteife bzw. steife Anfangswertaufgaben einschließlich Differential-Algebraischer Probleme enthalten. Wir verwenden das Matlab-Programm ode45, in demeingebettete Runge-Kutta-Verfahren der Ordnungen 4, 5 implementiert sind. Dieser furnichtsteife Anfangswertaufgaben geeignete Loser ist fur unser einfaches Anfangswertpro-blem ausreichend. Wir erinnern daran, daß wir bei der Nutzung der NAG-Bibliothek imAbschnitt 2 ebenfalls auf eingebettete Runge-Kutta-Verfahren der Ordnungen 4, 5 zuruck-gegriffen hatten. Die vollstandige Dokumentation zu ode45 findet man in der Hilfe zuMatlab.Wir benotigen ode45 wesentlich im function file fwert.m zur Funktionswertberechnung derzu minimierenden Funktion, vgl. fwert.m in c). Der verwendete Aufruf

[t1,x]=ode45(’reseitdgl’,[a b],x0)wird dort durch Kommentare teilweise erlautert. Der Parameter ’reseitdgl’ steht fur das

14

function file reseitdgl.m, in dem die rechte Seite der Differentialgleichung bereitgestellt wird.Das Ergebnis [t1,x] des Aufrufs ist eine dreispaltige Matrix. Diese enthalt zeilenweise die furdie Anfangswertaufgabe ermittelten diskreten Losungswerte (ti, xi1, xi2), wobei xi1 ≈ x1(ti),xi2 ≈ x2(ti). Der erste dieser Losungswerte ist der Anfangswert der Anfangswertaufgabe,und die Anzahl und damit die Zeilenzahl der Matrix ist wegen der Schrittweitensteuerungnicht von vornherein bekannt. Mit x(length(x),:) wird z.B. der letzte Zeilenvektor der Teil-matrix x von [t1,x] angesprochen.

c)”Matlab-Programm“ zur Losung der Optimierungsaufgabe:

Das”Matlab-Programm“ besteht aus drei M-Files, dem script file reaktor.m (Hauptpro-

gramm) und den beiden function files fwert.m und reseitdgl.m (Unterprogramme) furden Funktionswert der zu minimierenden Funktion und die rechte Seite der Differenti-algleichung. Alle drei Files sind jeweils unter ihrem Namen gespeichert. Wir geben jetztden vollstandigen Inhalt der drei Files an, versehen mit ausfuhrlichen Kommentaren. ImAnschluß daran gehen wir noch einmal auf wesentliche Gesichtspunkte ein:

% *****************************************************************************

% Scriptfile reaktor.m (Optimierung einer chemischen Reaktion) fuer Matlab 6.5

% *****************************************************************************

%

global n k10 k20 t % globale Parameter (in reaktor.m und im

global k1 k2 % function file fwert.m werden alle angegebenen

% Parameter gebraucht; im function

% file reseitdgl.m werden nur k1, k2 benoetigt)

n=input(’n= ’); % Eingabe n

h=1/n;

t=[0:h:1]; % Gitter t[0], ... , t[n] fuer Steuerfunktion

k10=480*0.892e9; % Konstanten

k20=480*0.768e16;

u0=ones(n,1).*350; % Startvektor fuer Minimierer

options=optimset(’TolFun’,1e-6,’TolX’,1e-6,’LargeScale’,’off’,’Display’,’iter’);

% Abbruchschranken fuer Minimierer (voreingestellt

% ist 1e-4); LargeScale ausschalten: Nutzer muss

% keine Ableitungen bereitstellen; Display iter:

% Ausgabe der Iterierten

%

[uopt,fopt]=fminunc(’fwert’,u0,options)

% Minimumstelle uopt und zugehoeriger Funktionswert

%

% Beginn Grafik ---------------------------------------------------------------

figure(1); clf % clf (clear figure) loescht den Inhalt des

% aktuellen Grafikfensters - hier des Fensters

% zu figure(1). Sonst wird beim Neustart

% von reaktor.m die neue figure(1) ueber die

15

% alte geschrieben.

figure(2); clf

figure(1) % Bild 1: Loesung der Differentialgleichung

% zu uopt

axis([0,1,0,0.8]) % Achsen und Skalierung

hold on % Bild fuer weitere Eintragungen erhalten

x0=[0.53 0.43]’; % Anfangswert der Differentialgleichung

for j=1:n

k1=k10*exp(-9000/uopt(j)); % k1,k2 fuer Intervall [t(j),t(j+1)]=:I(j)

k2=k20*exp(-15000/uopt(j));

a=t(j); % [a,b]=I(j)

b=t(j+1);

[t1,x]=ode45(’reseitdgl’,[a b],x0);

% Loesung der Dgl. auf [a,b] zum Anfangswert x0

plot(t1,x(:,1),’yo’,t1,x(:,2),’r’)

% Plotten der Loesung auf [a,b]

x0=x(length(x),:)’; % Endwert der Loesung auf [a,b] ist Anfangswert

% fuer naechstes Teilintervall

end

xlabel(’Zeit’)

title(’Konzentrationen’)

%

figure(2) % Bild 2: uopt=uopt(t) optimale Temperatur

axis([0,1,335,342])

hold on

for j=1:n

t2=[t(j) t(j+1)];

u2=[uopt(j) uopt(j)];

plot(t2,u2) % Plotten von uopt(t) auf Intervall I(j)

end

xlabel(’Zeit’)

title(’Temperatur’)

%

print Konzentrationen.eps -deps -f1

print Temperatur.eps -deps -f2

% Ausgabe der Bilder als eps-Files

% Option deps: Encapsulated PostScript

% Option depsc: Encapsulated Color PostScript

% Ende script file reaktor.m

% *****************************************************************************

% function file fwert.m (zu minimierende Funktion) fuer Matlab 6.5

% *****************************************************************************

16

function f=fwert(u)

global n k10 k20 t

global k1 k2

x0=[0.53 0.43]’; % Anfangswert der Differentialgleichung

for j=1:n

k1=k10*exp(-9e3/u(j)); % k1,k2 fuer Intervall [t(j),t(j+1)]=:I(j)

k2=k20*exp(-15e3/u(j));

a=t(j); % [a,b]=I(j)

b=t(j+1);

[t1,x]=ode45(’reseitdgl’,[a b],x0);

% Loesung der Dgl. auf [a,b] zum Anfangsw. x0

% bei voreingestellter Genauigkeit

% options=odeset(’AbsTol’,1e-7)

% [t1,x]=ode45(’reseitdgl’,[a b],x0,options);

% Loesung ... mit z.B. absoluter Genauigkeit 1e-7

x0=x(length(x),:)’; % Endwert der Loesung auf [a,b] ist Anfangswert

% fuer naechstes Teilintervall

end

f=-x0(2); % Funktionswert

% Ende function file fwert.m

% *****************************************************************************

% function file reseitdgl.m (rechte Seite der Differentialgleichung) fuer

% Matlab 6.5

% *****************************************************************************

%

function f=reseitdgl(t,x)

global k1 k2

f=zeros(2,1);

f(1)=-k1*x(1);

f(2)=k1*x(1)-k2*x(2);

% Ende function file reseitdgl.m

Im scriptfile reaktor.m und den function files fwert.m sowie reseitdgl.m werden globale Va-riablen verwendet. Auf eine globale Variable kann von jedem M-File, in dem die Variableglobal vereinbart ist, zugegriffen werden. Dies entspricht der Arbeit mit den COMMON-Bereichen bei der Fortran-Implementierung in Teil 2).In der zweiten Halfte von reaktor.m wird die ermittelte Losung der Optimierungsaufgabegraphisch veranschaulicht. Es werden zwei Bilder erzeugt. Im Bild 1 werden die der ermit-telten Minimumstelle uopt bzw. die der zugehorigen Treppenfunktion T (.;uopt) entspre-chenden Losungskomponenten x1(t) und x2(t) der Anfangswertaufgabe dargestellt. Dabeiwerden die Losungskomponenten x1(t) und x2(t) entsprechend dem im Abschnitt 1 be-schriebenen Anstuckeln der Losung zuerst im Intervall [t(1), t(2)] erzeugt und gezeichnet,

17

anschließend im Intervall [t(2), t(3)] usw. (vgl. j-Schleife fur Bild 1). Das Kommando holdon sichert, daß jeweils in das alte Bild gezeichnet wird. Andernfalls wurde jedem plot-Kommando ein neues Bild zugeordnet. Bei der Losung der Anfangswertaufgabe auf demjeweiligen Teilintervall wird wie bereits im Teil 2 wegen der Einfachheit des Problemsdavon ausgegangen, daß der Anfangswertloser erfolgreich arbeitet. Dies mußte eigentlichuberpruft werden. Bild 2 enthalt die Treppenfunktion T (.;uopt). Diese wird ebenfalls in-tervallweise gezeichnet. Bei der Abarbeitung von reaktor.m (vgl. d) erscheinen die beidenBilder in separaten Fenstern auf dem Bildschirm. Am Ende des Graphikteils von reaktor.mwerden die beiden Bilder als Encapsulated PostScript-Files Konzentrationen.eps und Tem-peratur.eps ausgegeben, so daß sie spater verfugbar sind. Selbstverstandlich konnten diebeiden Bilder weiter vervollkommnet werden, z.B. fehlen im Bild 1 die Funktionsbezeichnerx1(t), x2(t).Im function file fwert.m zur Funktionswertberechnung der zu minimierenden Funktion wirddas Anstuckeln der Losung der Anfangswertaufgabe mit der j-Schleife um den Aufruf desAnfangswertlosers realisiert. Auch hier wird stillschweigend davon ausgegangen, daß derLoser bei jedem Aufruf erfolgreich arbeitet.

d) Ausfuhren von reaktor.m: Der script file reaktor.m und die function files fwert.mund reseitdgl.m stehen im aktuellen Direktorie. Die Abarbeitung von reaktor.m erfolgtdurch das Kommando reaktor. Nach der Eingabe von n, als Beispiel wurde vier gewahlt,erscheinen fast unmittelbar die Ergebnisse:

>> reaktor

n= 4

Directional

Iteration Func-count f(x) Step-size derivative

1 4 -0.445182 1 -0.000379

2 13 -0.678808 1363.72 1.57e-05

3 22 -0.679468 1.52255 -2.09e-08

Optimization terminated successfully:

Current search direction is a descent direction, and magnitude of

directional derivative in search direction less than 2*options.TolFun

uopt =

339.1433

337.7757

336.9485

336.4483

fopt =

-0.6795

>>

18

Die Bilder 1 und 2 werden auf dem Bildschirm in separaten Fenstern ausgegeben. Im Bild1

”Konzentrationen“ ist die fallende Funktion x1(t), d.h. die Konzentration des Stoffes A1.

Die wachsende Funktion ist x2(t), d.h. die Konzentration des Stoffes A2:

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Zeit

Konzentrationen

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1335

336

337

338

339

340

341

342

Zeit

Temperatur

19

4 Bearbeitung der Optimierungsaufgabe mit Mathe-

matica

Eine kurze Einfuhrung in Mathematica findet man z.B. in [8]. Die Arbeit mit Mathematicawird durch die integrierte umfangreiche Hilfe und Dokumentation wesentlich unterstutzt. Inden folgenden Unterpunkten a) und b) gehen wir zunachst auf den Minimierer und den An-fangswertloser in Mathematica ein. Der anschließende Punkt c) enthalt das Mathematica-Programm (Notebook) zur Losung der Optimierungsaufgabe fur die Mathematica-Version5. Dieses Programm lauft nicht in der Vorgangerversion 4.2, vgl. Aufruf des Minimierersim folgenden Unterpunkt a). Ein fur die Version 4.2 modifiziertes Programm geben wir ind) an. Der Verfasser dankt Herrn Dr. Hans-Peter Scheffler, Fachrichtung Mathematik derTU Dresden, fur Hinweise zur Arbeit mit Mathematica, besonders fur die automatischeErzeugung des Minimierer-Aufrufs in der Version 4.2.

a) Minimierer: Mit dem Mathematica-Kommando FindMinimum konnen Funktionenvon mehreren Variablen ohne Nebenbedingungen numerisch minimiert werden.Wir beziehen uns zunachst auf die Mathematica-Version 5. Der von uns verwendete Aufruflautet

FindMinimum[f [u], {u, stavek}].Dabei ist f die zu minimierende Funktion, u der Variablenvektor und stavek der Startvek-tor fur das verwendete Minimierungsverfahren. Die Dimension des Variablenvektors u istdurch die Dimension von stavek festgelegt. FindMinimum hat weitere optionale Parame-ter, z.B. fur das Minimierungsverfahren, den Gradienten der zu minimierenden Funktionund Genauigkeitsparameter. Beim obigen Aufruf wird mit Voreinstellungen gearbeitet. Eswird z.B. ein Quasi-Newton-Verfahren verwendet, fur das der Nutzer keinen Gradientenbereitstellen muß. Bei erfolgreicher Arbeit liefert FindMinimum die numerische Losung inder Form {fmin, {u→ umin}}.FindMinimum ist in der Version 5 gegenuber der Version 4.2 in verschiedener Hinsicht ver-bessert worden. So wurden z.B. die einsetzbaren Minimierungsverfahren aktualisiert. Des-weiteren wurde der Aufruf vereinfacht. Obiger Aufruf ist in der Version 4.2 nicht moglich.In der Version 4.2 lautet der entsprechende Aufruf z.B. im Fall der Dimension n = 3

FindMinimum[f [{u1, u2, u3}], {u1, {st11, st12}}, {u2, {st21, st22}}, {u3, {st31, st32}}].Dabei ist unwesentlich, daß hier fur ableitungsfreie Rechnung ein zweiter Startwert verlangtwird. Im Aufruf sind z.B. st11 und st12 die beiden Startwerte fur die Variable u1. Wesent-lich ist, daß sich der Aufruf mit der Dimension n, d.h. der Anzahl der Variablen, andert.Wenn wie bei unserer Aufgabe die Dimension n nicht von vornherein fest ist, sondern alsEingangsparameter vom Nutzer gewahlt wird, ist dies wenig nutzerfreundlich. Man kannallerdings den Aufruf automatisch erzeugen, z.B. wie in d) realisiert.

b) Anfangswertloser: Zur numerischen Losung von Anfangswertaufgaben bei gewohnli-chen Differentialgleichungssystemen dient das Mathematica-Kommando NDSolve. Mit demKommando konnen auch ausgewahlte Probleme bei partiellen Differentialgleichungen undAlgebrodifferentialgleichungen gelost werden, was aber fur uns keine Rolle spielt. NDSolveist in der Mathematica-Version 5 gegenuber der Version 4.2 deutlich verbessert worden,

20

was sich jedoch bei der Anwendung auf unsere einfache Anfangswertaufgabe nicht auswirkt.Wir verwenden den Standardaufruf von NDSolve, der in beiden Versionen derselbe ist. Die-ser ist aus den Mathematica-Programmen in c) und d) ersichtlich, der Ruf von NDSolveerfolgt innerhalb der Funktion awa. Das Ergebnis von NDSolve - vgl. Probe fur Aufruf vonawa in c) - ist ein sogenanntes InterpolatingFunction object, das im Prinzip eine Nahe-rungsfunktion fur die Losung der Anfangswertaufgabe darstellt. Diese Naherungsfunktionwird durch Anwendung eines numerischen Diskretisierungsverfahrens, das Naherungswer-te fur die Losung der Anfangswertaufgabe an gewissen Stellen liefert, und anschließendegeeignete Interpolation erzeugt. Hier besteht ein Unterschied zu den Anfangswertlosernder NAG-Bibliothek und von Matlab, wo von den durch das Diskretisierungsverfahrenerhaltenen diskreten Naherungswerten nicht automatisch zu einer Naherungsfunktion furdie Losung der Anfangswertaufgabe auf dem gesamten Intervall, auf dem die Losung ge-sucht ist, ubergegangen wird. Dies bringt sicherlich Rechenzeitnachteile mit sich. Fur einComputeralgebra-System hat dies aber auch Vorteile. Beim Standardaufruf von NDSol-ve werden Voreinstellungen fur eine Reihe optionaler Parameter benutzt. So ist z.B. dieVerwendung von Adams-Verfahren fur nichtsteife Probleme und die Verwendung der BDF-Formeln (Gear) fur steife Probleme voreingestellt. Die Umschaltung nichtsteif/steif erfolgtautomatisch. Wegen der Einfachheit unserer Anfangswertaufgabe durfte bei uns nur Adamszum Einsatz kommen.Es sei erwahnt, daß der Quelltext von NDSolve 1400 Seiten lang ist, was wegen des breitenLeistungsspektrums nicht verwundert.

c) Mathematica-Programm (Notebook) fur Version 5 und Abarbeitung: DasNotebook fur die Version 5 ist auf den folgenden drei Seiten wiedergegeben. Es ist imFile Version5.nb gespeichert. Das Notebook enthalt in der vorliegenden Form nicht nurAnweisungen (Inhalt der In[ ]-Zellen), sondern auch die nach Abarbeitung des Notebookserzielten Ergebnisse. Ferner sind von In[4] bis In[6] und von In[8] bis In[13] Probeaufrufeund Anweisungen aufgenommen worden, die nur dem besseren Verstandnis dienen sollen.Der wesentliche Inhalt beginnt mit der Festlegung von Konstanten in den Zellen 1 und2. Diese Konstanten sind im gesamten Notebook verfugbar, also global. In Mathemati-ca sind alle Großen global bis auf die innerhalb einer Mathematica-Funktion vereinbar-ten Großen. So sind z.B. k1 und k2 in der in Zelle 3 definierten Mathematica-Funktionawa nur innerhalb von awa gultig, d.h. lokal. Die Mathematica-Funktion awa[a Real,x1aReal,x2a Real,b Real,u Real] liefert eine numerisch ermittelte Naherungsfunktion fur dieLosung der Anfangswertaufgabe x1 = −k1x1, x2 = k1x1 − k2x2, a ≤ t ≤ b, x1(a) =x1a, x2(a) = x2a, wobei k1 = k10e−9000/u, k2 = k20e−15000/u, vgl. auch b). Dabei wirdwie bereits bei der Arbeit mit der NAG-Bibliothek und mit Matlab stillschweigend an-genommen, daß es wegen der Einfachheit der Differentialgleichung keine Problem bei dernumerischen Losung gibt. Mit awa ist das im Abschnitt 1 beschriebene Anstuckeln derLosung fur die Anfangswertaufgabe (3) zu gegebener Treppenfunktion T (.;u) sehr einfach.In der in der Zelle 7 bereitgestellten Mathematica-Funktion f[u /;VectorQ[u,NumberQ]]zur Funktionswertberechnung der zu minimierenden Funktion wird in der Do-Schleife dieLosung angestuckelt. Die Mathematica-Funktion FM[stavek List,dim Integer] in der Zel-

21

le 14 ruft im Fall der Dimension dim=n den Minimierer bei Verwendung des Startvektorsstavek und stellt die Ergebnisse graphisch dar. Im Graphikteil von FM verhindert die Plot-Option DisplayFunction→Identity die Ausgabe der Graphik auf dem Bildschirm. Dagegenbewirkt DisplayFunction→$DisplayFunction, dies entspricht der Voreinstellung fur den op-tionalen Parameter DisplayFunction, Ausgabe auf dem Bildschirm. Mit dem KommandoShow werden zwei Graphiken zu einer vereinigt. Beim Aufruf des Minimierers wird wiederstillschweigend erfolgreiche Arbeit vorausgesetzt, d.h., daß eine

”Losung“ gefunden wird.

Beim Aufruf von FM in der Zelle 15 und bei der in der Zelle 13 vorangegangenen Probe istder Abbruch des Minimierers regular, die Textausgabe grenzt lediglich die Abbruchursacheein.Die Abarbeitung des Notebooks kann manuell zellenweise erfolgen oder in einem Zug uberKernel→Evaluation→Evaluate Notebook. Fur weitere Aufrufe von FM sind zusatzliche Zel-len analog zur Zelle 15 anzufugen, wobei der Startvektor stavek der Dimension anzupassenist.

22

Optimierung einer chem. Reaktion für Mathematica 5

� Konstanten

In[1]:= x10 = 0.53; x20 = 0.43

Out[1]= 0.43

In[2]:= k10 = 480*0.892*10^9; k20 = 480*0.768* 10^16

Out[2]= 3.6864´1018

� Funktion awa zur numerischen Lösung des Anfangswertproblems x 1 = −k1 x1 , x 2 = k1 x1 - k2 x2 , actcb, x1 (a)= x1 a , x2 (a) = x2 a ,wobei k1 = k10 e-9000�u, k2 = k20 e-15000�u und a, x1 a, x2 a , b sowie u Parameter .

In[3]:= awa@a_Real, x1a_Real, x2a_Real, b_Real, u_RealD :=Module@8k1 = k10* Exp@-9000�uD, k2 = k20*Exp@-15000�uD<,NDSolve@8x1’@tD � -k1* x1@tD, x2’@tD � k1*x1@tD - k2*x2@tD,x1@aD � x1a, x2@aD == x2a<, 8x1, x2<, 8t, a, b<DD

� Probe für Aufruf von awa *************************************************************************

In[4]:= l1 = awa@0., x10, x20, 100., 250.DOut[4]= 88x1 ® InterpolatingFunction@880., 100.<<, <>D,

x2 ® InterpolatingFunction@880., 100.<<, <>D<<In[5]:= 88hilf1, hilf2<< = 8x1@100.D, x2@100.D< �. l1 H* Lösung an der Stelle 100 *LOut[5]= 880.524762, 0.435236<<In[6]:= hilf1

Out[6]= 0.524762

********************************************************************************* Ende Probe awa

� Funktion f zur Berechnung des Funktionswertes der zu minimierenden Funktion.Wegen der Ergänzung /; VectorQ[u,NumberQ] im Argument wird f [u] beim Aufruf nur unter der Bedingung ausgewertet,dass u ein Vektor aus Zahlen ist. Sonst wird f[u] zurückgegeben. Die Auswertung von f[u] wäre dann auch gar nichtmöglich, weil dies der Anfangswertlöser awa nicht zuläßt. Die Ergänzung hat Bedeutung beim Aufruf von FindMinimum.Ohne die Ergänzung funktioniert der Aufruf nicht.

In[7]:= f@u_ �; VectorQ@u, NumberQDD := Module@8x1e = x10, x2e = x20, i<,Do@Hl1 = awa@tvek@@iDD, x1e, x2e, tvek@@i + 1DD, u@@iDDD;88x1e, x2e<< = 8x1@tvek@@i + 1DDD, x2@tvek@@i + 1DDD< �. l1L, 8i, 1, n<D;-x2e D

� Probe für Aufruf von f und FindMinimum ***********************************************************

In[8]:= n = 4; h = 1�n; tvek = N@Table@Hi - 1L h, 8i, 1, n + 1<DDOut[8]= 80., 0.25, 0.5, 0.75, 1.<

23

In[9]:= stavek = 8350., 350., 350., 350.<Out[9]= 8350., 350., 350., 350.<In[10]:= f@stavekDOut[10]= -0.445182

In[11]:= f@8para, 350., 350., 350.<DH* Aufruf wird zurückgegeben, weil Argument kein Vektor aus Zahlen *LOut[11]= f@8para, 350., 350., 350.<DIn[12]:= sol = FindMinimum@f@uD, 8u, stavek<D

FindMinimum::lstol : The line search decreased the step size to within tolerance specified by AccuracyGoal and

PrecisionGoal but was unable to find a sufficient decrease in the function. You may needmore than MachinePrecision digits of working precision to meet these tolerances. Mehr¼

Out[12]= 8-0.679488, 8u ® 8339.487, 337.811, 336.822, 336.162<<<In[13]:= First@sol@@2DDD@@2DD H* Minimumstelle *LOut[13]= 8339.487, 337.811, 336.822, 336.162<******************************************************************** Ende Probe f und FindMinimum

� Funktion für Aufruf von FindMinimum und grafische Darstellung der ermittelten Steuerfunktion u(t) und der zugehörigenKonzentrationen x1 (t) sowie x2 (t) im Fall der Dimension dim mit Startvektor stavek für Minimierer.

In[14]:= FM@stavek_List, dim_IntegerD :=Module@8h, g1, g2, aw<,n = dim; h = 1�n; tvek = N@Table@Hi - 1L h, 8i, 1, n + 1<DD;sol = FindMinimum@f@uD, 8u, stavek<D; Print@solD; v = First@sol@@2DDD@@2DD;H* Bild der Steuerfunktion *Lg1 = Plot@v@@1DD, 8t, tvek@@1DD, tvek@@2DD<, DisplayFunction ® IdentityD;Do@H g2 = Plot@v@@kDD, 8t, tvek@@kDD, tvek@@k + 1DD<, DisplayFunction ® IdentityD;

g1 = Show@g1, g2DL, 8k, 2, n<D;Show@g1, AxesLabel ® 8t, u<, DisplayFunction ® $DisplayFunctionD;H*Bild für Komponenten der Zustandsfunktion*Ll1 = awa@tvek@@1DD, x10, x20, tvek@@2DD, v@@1DDD;g1 = Plot@Evaluate@8x1@tD, x2@tD< �. l1D,8t, tvek@@1DD, tvek@@2DD<, DisplayFunction ® IdentityD;Do@Haw = First@8x1@tvek@@kDDD, x2@tvek@@kDDD< �. l1D;l1 = awa@tvek@@kDD, aw@@1DD, aw@@2DD, tvek@@k + 1DD, v@@kDDD;g2 = Plot@Evaluate@8x1@tD, x2@tD< �. l1D,8t, tvek@@kDD, tvek@@k + 1DD<, DisplayFunction ® IdentityD;g1 = Show@g1, g2DL, 8k, 2, n<D;

Show@8g1,Graphics@8Text@"x1HtL", 80.8, 0.32<, 80, 1<D, Text@"x2HtL", 80.8, 0.6<, 80, 1<D<D<,AxesLabel ® 8t, ""<, DisplayFunction ® $DisplayFunctionDD

� Aufruf für Dimension 4 mit dem Startvektor stavek

24

In[15]:= FM@stavek, 4DFindMinimum::lstol : The line search decreased the step size to within tolerance specified by AccuracyGoal and

PrecisionGoal but was unable to find a sufficient decrease in the function. You may needmore than MachinePrecision digits of working precision to meet these tolerances. Mehr¼8-0.679488, 8u ® 8339.487, 337.811, 336.822, 336.162<<<

0.2 0.4 0.6 0.8 1t

336.5

337.5

338

338.5

339

339.5u

0.2 0.4 0.6 0.8 1t

0.3

0.4

0.5

0.6

x1HtLx2HtL

Out[15]= � Graphics �

25

d) Mathematica-Programm (Notebook) fur Version 4.2: Die folgenden dreiSeiten enthalten das Notebook fur die Version 4.2. Dieses unterscheidet sich vom Note-book fur die Version 5 aus c) im wesentlichen im Aufruf des Minimierers, vgl. a). Biszur Mathematica-Funktion f[u /;VectorQ[u,NumberQ]] zur Funktionswertberechnung derzu minimierenden Funktion sind beide Notebooks

”identisch“, im Notebook fur die Versi-

on 4.2 wird in diesem Teil lediglich auf Kommentare, Ergebnisse und Proben verzichtet.Die Konstruktion des Aufrufs fur den Minimierer erfolgt mit der Mathematica-FunktioncallFindMinimum[f Symbol,stvek0 List,stvek1 List] in der Zelle 5. Die Arbeitsweise voncallFindMinimum wird schrittweise in den Zellen 6 bis 15 erlautert. In der MathematicaFunktion FM in der Zelle 16 wird dieser Aufruf genutzt. FM unterscheidet sich von derMathematica-Funktion FM in der Version 5 nur in diesem Punkt und daß jetzt zwei Start-werte erforderlich sind.

26

Optimierung einer chem. Reaktion für Mathematica 4.2

In[1]:= x10 � 0.53; x20 � 0.43;

In[2]:= k10 � 480 � 0.892� 10^9; k20 � 480 � 0.768 � 10^16;

In[3]:= awa�a_Real, x1a_Real, x2a_Real, b_Real, u_Real� :�Module��k1 � k10 � Exp��9000 � u�, k2 � k20 � Exp��15000� u��,NDSolve��x1'�t� � �k1 � x1�t�, x2'�t� � k1 � x1�t� � k2 � x2�t�,x1�a� � x1a, x2�a� �� x2a�, �x1, x2�, �t, a, b���

In[4]:= f�u_ �; VectorQ�u, NumberQ�� :� Module��x1e � x10, x2e � x20, i�,Do��l1 � awa�tvek��i��, x1e, x2e, tvek��i � 1��, u��i���;��x1e, x2e�� � �x1�tvek��i � 1���, x2�tvek��i � 1���� �. l1�, �i, 1, n��;

�x2e �

� Konstruktion des Aufrufs für Minimierer und Aufruf. Die Argumente stvek0 und stvek1 sind die Startvektoren fürden Minimierer.� Das Kommando Remove["u*"] löscht alle Variablen, die mit u beginnen. Die Indizierung der Komponenten desim darauf folgenden Kommando var=Table[Unique["u"],{k,1,n}] erzeugten Variablenvektors beginnt dann immer mit1. Siehe auch die Zellen 13 bis 15.�Die restlichen Kommandos werden anschließend erläutert.

In[5]:= callFindMinimum�f_Symbol, stvek0_List, stvek1_List� :�Module��steuervek, start�,Remove�"v�"�; steuervek � Table�Unique�"v"�, �k, 1, n��;

start � Table��steuervek��k��, �stvek0��k��, stvek1��k����, �k, 1, n��;sol � Apply�FindMinimum, Prepend�start, Apply�f, �steuervek����

�

� Erläuterung der Kommandos in callFindMinimum und Probe für Aufruf von callFindMinimum im Fall n=4

In[6]:= n � 4; h � 1 � n; tvek � N�Table��i � 1� h, �i, 1, n � 1���

Out[6]= �0., 0.25, 0.5, 0.75, 1.�

In[7]:= steuervek � Table�Unique�"v"�, �k, 1, n��

Out[7]= �v1, v2, v3, v4�

In[8]:= Table�Unique�"v"�, �k, 1, n��

Out[8]= �v5, v6, v7, v8�

In[9]:= Remove�"v�"�; steuervek � Table�Unique�"v"�, �k, 1, n��

Out[9]= �v1, v2, v3, v4�

In[10]:= stvek1 � �350., 350., 350., 350.�; stvek0 � �338., 338., 338., 338.�;

In[11]:= start � Table��steuervek��k��, �stvek0��k��, stvek1��k����, �k, 1, n��

Out[11]= ��v1, �338., 350.��, �v2, �338., 350.��, �v3, �338., 350.��, �v4, �338., 350.���

27

In[12]:= Apply�f, �steuervek���� Bei Apply�f,�steuervek�� wird zunächst f��v1,v2,v3,v4�� gebildet und

anschließend aufgerufen. Da das Argument kein Vektor aus Zahlen ist,wird der Aufruf nicht ausgewertet sondern zurückgegeben. ��

Out[12]= f��v1, v2, v3, v4��

In[13]:= Prepend�start, Apply�f, �steuervek���

Out[13]= �f��v1, v2, v3, v4��, �v1, �338., 350.��,�v2, �338., 350.��, �v3, �338., 350.��, �v4, �338., 350.���

In[14]:= Apply�FindMinimum, Prepend�start, Apply�f, �steuervek����

Out[14]= ��0.679489, �v1 � 339.486, v2 � 337.806, v3 � 336.83, v4 � 336.157��

In[15]:= callFindMinimum�f, stvek0, stvek1� �� Probe für Aufruf von callFindMinimum ��

Out[15]= ��0.679489, �v1 � 339.486, v2 � 337.806, v3 � 336.83, v4 � 336.157��

*********************************************************** Ende Erläuterung/Probe callFindMinimum

� Funktion für Aufruf von callFindMinimum und grafische Darstellung der ermittelten Steuerfunktion u(t) und derzugehörigen Konzentrationen x1 (t) sowie x2 (t) im Fall der Dimension dim mit Startvektoren stvek0 und stvek1 fürMinimierer.

In[16]:= FM�stvek0_List, stvek1_List, dim_Integer� :�Module��h, g1, g2, aw�,n � dim; h � 1 � n; tvek � N�Table��i � 1��h, �i, 1, n � 1���;callFindMinimum�f, stvek0, stvek1�; Print�sol�;

�� Bild der Steuerfunktion ��

g1 � Plot�sol��2����1����2��,�t, tvek��1��, tvek��2���, DisplayFunction � Identity�;

Do�� g2 � Plot�sol��2����k����2��, �t, tvek��k��, tvek��k � 1���,DisplayFunction � Identity�;g1 � Show�g1, g2��, �k, 2, n��;

Show�g1, AxesLabel � �t, u�, DisplayFunction � $DisplayFunction�;

��Bild für Komponenten der Zustandsfunktion��l1 � awa�tvek��1��, x10, x20, tvek��2��, sol��2����1����2���;g1 � Plot�Evaluate��x1�t�, x2�t�� �. l1�,�t, tvek��1��, tvek��2���, DisplayFunction � Identity�;

Do��aw � First��x1�tvek��k���, x2�tvek��k���� �. l1�;l1 � awa�tvek��k��, aw��1��, aw��2��, tvek��k � 1��, sol��2����k����2���;g2 � Plot�Evaluate��x1�t�, x2�t�� �. l1�,�t, tvek��k��, tvek��k � 1���, DisplayFunction � Identity�;

g1 � Show�g1, g2��, �k, 2, n��;Show��g1, Graphics�

�Text�"x1�t�", �0.8, 0.32�, �0, 1��, Text�"x2�t�", �0.8, 0.6�, �0, 1�����,AxesLabel � �t, ""�, DisplayFunction � $DisplayFunction�

�

28

� Aufruf für Dimension 4 mit den Startvektoren stvek0 und stvek1

In[17]:= FM�stvek0, stvek1, 4�

��0.679489, �v1 � 339.486, v2 � 337.806, v3 � 336.83, v4 � 336.157��

0.2 0.4 0.6 0.8 1t

336.5

337.5

338

338.5

339

339.5

u

0.2 0.4 0.6 0.8 1t

0.3

0.4

0.5

0.6

x1�t�

x2�t�

Out[17]= ��Graphics��

29

5 Rechenzeiten

a) Art der Zeitmessung:

NAG: Verwendung der Fortran-Standard-Subroutine CPU TIME. Sei im Hauptprogrammreaktor z.B. zeit REAL-Variable. Dann liefert der Aufruf CPU TIME(zeit) am Programm-ende die Rechenzeit zeit.

Matlab: Verwendung des Kommandos cputime. Im Skript reaktor.m werde am Anfangdie Anweisung zeit1=cputime; und am Ende die Anweisung zeit2=cputime; aufgenom-men. Dann ist die Differenz zeit2-zeit1 die Rechenzeit.

Mathematica: Verwendung des Kommandos Timing. Im Notebook fur die Version 5 er-gibt z.B. n = 5; Timing[FM [Table[350, {i, n}], n] die Rechenzeit im Fall n = 5.

b) Rechenzeiten und Vergleich:

NAG auf Workstation IBM RS/6000

n Zeit in Sekunden fmin

5 0.37 -0.67950710 1.23 -0.67953415 2.45 -0.67954020 4.21 -0.67954125 6.34 -0.67954230 8.94 -0.67954350 21.06 -0.67954375 52.50 -0.679544

100 94.65 -0.679544

Matlab 7 auf PC (2.4 GHz)

mit Minimierer fminunc mit Minimierer fminsearchn Zeit in Sekunden fmin n Zeit in Sekunden fmin

5 4.69 -0.679507 5 19.86 -0.67950710 14.34 -0.679534 10 145.68 -0.67953415 30.49 -0.679539 15 425.92 -0.67953920 52.62 -0.679541 20 758.89 -0.67954125 80.63 -0.679542 25 1180.35 -0.67954230 116.04 -0.679543 30 1721.08 -0.67954250 314.53 -0.67954375 680.88 -0.679543

100 1045.57 -0.679543

30

Mathematica auf PC (2.4 GHz)

Version 5 Version 4.2n Zeit in Sekunden fmin n Zeit in Sekunden fmin

5 3.14 -0.679507 5 1.94 -0.67950910 12.38 -0.679534 10 4.84 -0.67953015 26.87 -0.679540 15 17.6 -0.67953420 35.30 -0.679541 20 30.02 -0.67953525 68.30 -0.67954230 82.32 -0.67954350 413.29 -0.67954375 470.00 -0.679543

100 1408.43 -0.679544

20 40 60 80 100n

200

400

600

800

1000

1200

1400

tMathematica

Matlab:fminunc

NAG

Bild: Zeit t in Sekunden in Abhangigkeit von n

Die ermittelten Rechenzeiten zeigen (vgl. Bild und Tabellen), daß Matlab und Mathe-matica in etwa vergleichbar sind und NAG deutlich uberlegen ist. Bei dem Vergleich derRechenzeiten von Matlab und Mathematica einerseits und NAG andererseits ist zwar zubeachten, daß die Rechnungen mit Matlab und Mathematica auf einem PC (2.4 GHz) so-wie die mit NAG auf einer Workstation (IBM RS/6000, Baujahr 1996/97) durchgefuhrtwurden. Folglich ist ein absoluter Vergleich der Rechenzeiten schwer moglich. Entschei-dend ist aber, daß die Rechenzeit bei NAG weniger schnell mit n wachst als bei Matlabund Mathematica.Naturlich ist dieses Ergebnis zu erwarten. Denn Matlab und Mathematica arbeiten als In-terpreter. Hier wird Kommando fur Kommando interpretiert und abgearbeitet. Dies fuhrtzu einem starken Anwachsen des Aufwandes mit wachsendem n. Dagegen erzeugt bei NAGder Compiler ein Maschinenprogramm und gemessen wird die Laufzeit dieses Maschinen-programmes.

31

6 Schlußfolgerungen

Die Losung des Steuerproblems ist mit allen drei Softwaresystemen, d.h. mit der NAG-Bibliothek, mit Matlab und Mathematica, prinzipiell moglich.Fur das Computeralgebra-System Mathematica ist dies durchaus bemerkenswert. Mit Maplekonnte das Steuerproblem nicht bearbeitet werden, weil dort das Kapitel Optimierung un-zureichend vertreten ist. Generell ist festzustellen, daß die NAG und auch Matlab imnumerischen Rechnen den Computeralgebra-Systemen uberlegen sind. Die implementier-ten Verfahren sind haufig aktueller, und es sind mehr Aufgabenklassen vertreten. Aller-dings versuchen die Computeralgebra-Systeme fortlaufend, den Numerikteil auszubauen.Das Mathematica-Notebook fur die Version 4.2, konkret der umstandliche Aufruf des Mi-nimierers, offenbart konzeptionelle Probleme beim Umgang mit Numeriksoftware. DiesesProblem ist in der Version 5 behoben. Der Aufwand fur die Aktualisierung von Notebooksder Version 4.2 fur die Version 5 ist jedoch nicht zu vernachlassigen wie die Notebooksfur das Steuerproblem zeigen. Dagegen war bei der NAG und bei Matlab der Ubergangvon Version zu Version problemlos, weil dort hochstens ein Verfahren durch ein anderes zuersetzen war und an der Konzeption nichts geandert wurde.Ein wesentlicher Vorteil der NAG gegenuber Matlab und Mathematica ist die deutlichgeringere Rechenzeit bei großer Dimension n. Die Rechenzeit fur das bei der NAG vomCompiler erzeugte Maschinenprogramm wachst mit n wesentlich langsamer als die Rechen-zeiten fur die interpretierende Abarbeitung der Skriptfiles bei Matlab und Mathematica.Die unterschiedliche Arbeitsweise der drei Softwaresysteme zieht ferner nach sich, daßdas mit der NAG erzeugte Maschinenprogramm auf dem entsprechenden Rechnertyp oh-ne die NAG-Bibliothek lauffahig ist. Dagegen muß zur Abarbeitung eines Matlab- bzw.Mathematica-Files stets Matlab bzw. Mathematica selbst vorhanden sein.Ein Pluspunkt fur Matlab und Mathematica ist die integrierte Graphiksoftware. Diese er-leichtert die Darstellung der Ergebnisse. Ebenso kann das in Mathematica und Matlabintegrierte symbolische Rechnen von Vorteil sein. Bei Matlab sind dazu allerdings mogli-cherweise die Symbolic Math Toolbox oder die Extended Symbolic Math Toolbox erforderlich,die nicht im Lieferumfang des Grundpaketes enthalten sind. Bei der NAG muß fur graphi-sche Darstellungen und zur Einbeziehung symbolischer Rechnungen das System gewechseltwerden.

32

Literatur

[1] W.Alt: Nichtlineare Optimierung. Eine Einfuhrung in Theorie, Verfahren und Anwen-dungen. Vieweg 2002.

[2] L.Cesari: Optimization - Theory and Applications. Springer, New York 1983.

[3] Yu.G.Evtushenko: Methoden zur Losung von Extremwertaufgaben und ihre Anwen-dung in Optimierungssystemen. (Metody resheniya ehkstremal´nykh zadach i ikhprimenenie v sistemakh optimizatsii. Reihe: Optimizatsiya i Issledovanie Operatsij.)Nauka, Moskva 1982. (In Russisch)

[4] L.T.Fan: The Continuous Maximum Principle. John Wiley & Sons, New York 1966.

[5] E.Hofer, R.Lunderstadt: Numerische Methoden der Optimierung. Oldenbourg,Munchen Wien 1975.

[6] H.J.Oberle, G.Opfer: Optimale Steuerung, Einfuhrung in die Theorie und numeri-sche Verfahren (Teil 2). Hamburger Beitrage zur Angewandten Mathematik, Reihe B,Bericht 13, 1989.

[7] L.S.Pontrjagin, V.G.Boltjanskij, R.V.Gamkrelidze, E.F.Miscenko: MathematischeTheorie optimaler Prozesse. Oldenbourg, Munchen 1967.

[8] W.Preuß, G.Wenisch (Herausgeber): Lehr- und Ubungsbuch Numerische Mathematik.Fachbuchverlag Leipzig im Carl Hanser Verlag 2001.

[9] H.-G.Roos, H.Schwetlick: Numerische Mathematik. Reihe Mathematik fur Ingenieureund Naturwissenschaftler. B.G.Teubner, Stuttgart/Leipzig 1999.

[10] J.Stoer, R.Bulirsch: Numerische Mathematik 2. Springer, funfte Auflage 2005.

[11] K.L.Teo, C.J.Goh, K.H.Wong: A unified computational approach to optimal controlproblems. Longman Scientific & Technical 1991.

33