Optimaliseringsmethoden en technieken in de technologie

DEC. 1979

ARCH lEE

Lab. y. Scheepsbouwkunde

lechnische Hogeschool

Technische Hogeschao Deift

Onderafdeling der Wiskunde

Vakgroep Toegepaste Analyse

Optimaliseringsmethoden

en - technieken

in de Technologie

door

Dr. C. de Wit

VakgroepToegepaste Analyse

Opticialiseringsmethoden en - technieken in de Technologie door Dr. C. de .Jit Collegedictaat a 165 A , 1977/78

-1-Inhoud.soDgave. _Biz.

O. Voorwoord en samenvatting. ₃

1 Introductie van het begrip "optimaliseren".

Voorbeelden0 ₇

1.1. Statisch en Dynamisch Optimaiiseren0 ₇

l.2 _{Voorbeelden van Statisch Optimaliseren0}

9 1.3e _{Enige wiskundige notaties, definities en}

conventies. li

2 _{Minimalisering van functies zonder beperkingen.} 13 2cl. Nodige en voidoende voorwaarden voor een vrij

locaal minimumpunt0 ₁₃

22

Enige practische voorbereidingsmaatregelen0 ₁₇

22l. Schaling van de variabelerì0 ₁₇

2.22

Numeriek djfferentjren0 ₁₈

2.3e Lijnminimaiisaties0 ₂₁

2.4e _{Minimaiiseringsmethoden, waarbj geen gebruik}

wordt gemaakt van gradinten0 ₂₃

24l De Directe Zoekmethod.e van Hooke en Jeeves0 23

242 Neider & Nead's Simplex Algoritme0 ₂₇

2.43. Kwad.ratisch convergente minimaliseringsmethode

van Powell0 ₃₀

2 5 Gradint-minimaliseringsmethoden0 32

25l. De anti-gradient of steilste dalingsmethode. ₃₂

252 De geconjugeerde gradintmethode0

₃₄

253 De methode van Newton.

₃₈

25.4

Quasi-Newton methoden0 ₄₁

2.5.5e Kort overzicht van het, op het

3,2.

Gelijkheidsbeperkingen,

47

3.30

Ongelijkheidsbeperkingen.

50

3.40

Straffuncties,

53

3.4.1. Inleiding.

53

3,4.2e Uitwendge straffuncties.

55

3,4,3. Voorbeeld van een SUMT-programiia,

57

3.4.4. Inwendige straffuncties.

60

305e

Verwerking van restricties door middel van

stijgende afknottingsniveau's.(Lit, 12)

62

3.5,1, Staha's COMET-methode,

62

3,5,2. Voorbeeld van een

.functie-rninimalisatiepro--gramma volgens de COMET-methode0

67

3,60

Minimalisering van functies onder lineaire

(on)geli,jkheidsbeperkingen, (Lit. i7)

70

4.

Vraagstukken.

72

4.1.

Minimalisering van fimcties zonder beperkingen.

72

4,2,

Minimalisering van functies met beperkingen.

76

O. Voorwoord en Samenvatting.

In dit college zullen diverse optimaliseringsmethoden en

-technieken worden besproken, die in de laatste decennia zijn

ontwikkeld. De nadruk zal daarbij vallen op het oplossen van die

technische problemen, waarvan bet mathematische model een

duide--l1ke niet-lineaire structur heeft. Voorts is bet betrekkelijk

geringe aantal variabelen - meestal minder dan tien en nooit

meer dan vijftig - min of meer karakteristiek voor de hier aan

de orde zLjnde technische optimaliseringsproblemen.

Indien de te optimaliseren functie voortkomt uit het

mathema--tisch modelleren van bet

n of andere

t.echnisch-fysische of

hemische gebeuren, dan mag worden verondersteld, dat de in het

geding zãjnde functies "voldoende glad" zijn, d..z. minstens

twee keer continu differentieerbaar naar al hun variabelen. Wie

dit een enigszins star uitgangspunt vindt, bedenke dat bij

veer--beeld functies als

1(x)

=

of

g(x) = sign(x) = x/jx(

glad te maken zijn door ze te vervangen door

2 2

f (x) =

(x

_+E)

en

g (x) = x/(x

_{+ E)}

met een voldoend kleine, maar wel positieve c

De, steeds regle, waarden van de te optimaliseren functies

.nnen in diverse vorrnen worden gereproduceerd. Vaak zln ze in

gesloten analytische vorm gegeven. Het is dan een betrekkeljk

eenvoudige zaak, in

nzelfde procedure naast de functiewaarden

ook de n-vector van partiale afgeleiden te construeren. Bj

der--gelijke functies ugt het betrekkelijk voor de hand, gebruik te

maken van optirnaliseringsmethoden, die veelvuldig met zowel

I'unctiewaarden als gradinten werken.

In toenemende mate echter komt het voor, dat de functiewaarden

alleen maar kunnen worden geproduceerd door middel van een

numerieke procedure. Nen denke h.v. aan de oplossingen van het

stelsel van niet-lineaire differentiaalvergeljkingen

i = 1(l)n

:

.(t) = f(x1(t),

. . _.

,x(t),a1,

. . _.

,a)

anuit een gegeven startpunt x(0)

De waarden x.(l)

,

i = l(l)n

,

zullen nu van de parametervector

Indien flu die f1 s continu differentieerbaar van a afhangen, dan

is deze F zeker twee keer continu di±'ferentieerbaar naar de

variabelen a1 t/m a. Deze partile afgeleiden zin echter niet

analytisch te produceren, evenmin als F(a) zeif.

'1el is het

moge1ijk, die partile afgeleiden te produceren via numerieke

diferentiatie van F. De hiervoor 'cenodigde waarden van

kunnen worden verkregen door nunerieke integratie van het

gege--ven stelsel van differentiaalvergelijkingen.

Wanneer men zo'n

F

wil minimaliseren, dan zal het duidelk zijn, dat men

veel

meer is genteresceerd in

optinaliseringsmethoden, die zo

weinig moge1jk func tiewaardeberekeningen vra;en.

:uflctiewaarden kunnen, behalve door een mathematisch model,

cok geleverd worden door registratie van het

n of andere

fysische, mechanische of electrische verschijnse1.

Dergel1jke

functiewaarden worden meestal dermate door ruis verstoord,

dat de verkregen beonstering niet differentieerbaar

is.

'Ien is

dan aangewezen op methoden, die deze differentieerbarheid

niet voorveonderstellen.

Het zal duideli,ik zijn, dat er, gelet op de hierboven

genoemde

roge1ikheden, niet zoiets bestaat als een "optimale"

optima--liseringsmethode. Men zal, gelet op ervaringen met bestaande

technieken, telkens opnieuw moeten schatten en proberen,

welke

methode er voor het onderhavige geval als de beste uitkomt.

Naast de keuze van de methode is er steeds een

keuzemoeiLjkheid

ten aanzien van het al dan niet gebruiken van eventuele

soft--ware. Voor het

nmalig oplossen van een niet al te groot

probleem is het vaak gemakkelijk en verleidelijk,

gebruik te maken

van voorgeprograimeerde optimaliseringsprocedures. Het voordeel

hiervan is natuur1ijk de betrekkeljk eenvoudige programmering year

de gebruiker. Hiertegenover echter staat het nadeel, dat men

bij een eventuele foutmelding niet kan nagean, waar

deze fout in

de Thlack-box-procedur&' precies optreedt. Men

kan hier in feite

alleen uitkomen, door van zo'n procedure een

listing op te vragen,

die in een eigen dataset overnernen en daarmee

de zaak nag eens

5

optreden, kan het zeifs beter zijn, een eigen, op het probleern

gerichte versie te maken van zo'n optimaliseringsmethode.

Tegen--over de moeite van het programmeerwerk staat dan het gemak van

de veel grotere flexibiliteit van zo'n programma.

Gelet op deze zaken zullen de meest voorkomende

optimaliserings--procedures in dit collegedictaat zo goed mogeUjk worden

beschre--ven. Tevens zal worden vermeld, op welke wijze deze methoden in

1gol aanroepbaar zijn, welke voorzorgsmaatregelen men moet

treffen en welke moeiLjkheden men kan verwachten bij het gebruik

van deze methoden.

Bu

het minirnaliseren van een objectfunctie treden vrjwel steeds

beperkingen op ten aanzien van de variabelen, waar die objectfunctie

van af hangt. Grootheden als de diameter van een scheepsschroef

of het verwarmend opperviak van een stoomketel hebben zekere,

door technologische overwegingen bepaalde, grenzen. Ook kan het

voorkomen, dat eon functie van meerdere variabelen niet boyen of

onder bepaalde grenzen mag komen. Hen denke bi

voorbeeld aan het

feit, dat een buigspanning van een constructie van een bepaald

materiaal een zekere, voor dat materiaal vastgestelde, maximale

waarde niet mag overschri,jden.

Deze onge1jkheidsbeperkingen (Unequality constraints) kunnen

worden geformuleerd als

g(x)O

,

j = l(l)m

.

(5.1)

Verder blijkt ook vaak, dat er bepaalde gelijkheidsrelaties tussen

de variabelen bestaan. Deze relaties kunen o.m. voortkomen

uit de eis, dat er in een bepaalde stationaire toestand

even--wicht moet zijn in mechanische, thermische of electrische zin.

Deze gelijkheidsrestrìcties (Equality constraints) kunnen worden

geformuleerd als

hk(x) = O

,

k = l(l)p <n

_.

_(5.2)

}iet totale probleem

T'Minimaliseer f(x), waarbj x moet voldoen aan

O _,

j = 1(l)m

hk(x) = O

,

k = 1(l)p ,"

kan op diverse wijzen worden opgelost. In dit college worden

alleen die oplosingsmethoden behandeld, die voldoende zijn

In- en uitwend.ige straffunctiemethoden,

Iitwendige afknottingsmethode, (iii.) Gradint-projectiemethode.

Als sluits-tuk van dit dictaat zijn een aantal vraa'stukken

opgenomen, die met behuip van de diverse behandelde

technieken en gedeeltelijk oak via de voorgeprogrammeerde procedures van het T0H0-Rekencentrum kunnen worden

-7

1.

_{Introductie van het begrip tOptimaliseren". Enige}

concrete voorbeelden.

1.1. 3tatisch en Dynamisch Optirnaliseren.

Vrj algemeen kan men in de technologie een

optimalise--ringsprobleem stellen, als een technisch

constructie-of besturingsprobleem meerdere oplossingen toelaat.

Bij een constructieprobleem kan men vaak een aantal

variabelen x1 t/m x

binnen zekere grenzen vrj kiezen.

Kiest men nu een reelwaardige, van onderen begrenade

functie f van die variabelen - voor zo'n f wordt meestal

een kosten- of verliesfunctie genornen - dan kan men als

optimaliseringsprobleem stellen:

"Kies de variabelen x1 t/m x

binnen de vanuit

de techniek volgende grenzen zo, dat f(x1,..,x)

rninimaal wordt."

bien spreekt in dit geval van "Statisch Optirnaliseren".

Bij besturingsproblemen ugt de zaak duidelijk anders.

en heeft in zo'n geval te maken met een aantal

fase-of toestandsvariabelen (State variables), dIe met de

tjd veranderen

_{: x1(t) t/m x(t) .}

_{Deze variabelen}

veranderen met de tjd als gevoig van

de waarde, die ze op elk moment hebben en t.g.v.

de waarde van

,

op het systeem van buitenaf werkende,

besturingsgrootheden (control variables) u1(t) t/m um(t)

Voorts kan dit dynamische gedrag nog expliciet van de

tijd afhangen. Er geldt dus

Jt) = f(x(t), u(t), t)

,

i = l(l)n

_(7.1)

x1(t)

_u1(t)

met x(t) =

:

en

u(t) =

-

xt)

_-

_u(t)

Eij zo'n dynamisch systeem is het altijd mogelijk, de

fase--vector x(t) vanaf een zekere begintoestand X0

O

tijd

to op meer dan

n manier over te brengen naar de

toestand

X o

tijd t1. Dit kan door goede keuze van de

Het is in zo'n geval zinvol, een naar beneden begrensde verlies- of kostenaangroeifunctie L te definiren, werkend op x(t), op u(t) en eventueel op t zeif:

L(x(t),u(t),t) L0

Men beschouwt nu binnen de klasse van toelaatbare stuur--functies u de daarb&j behorende waarde van de integraal

ti

=

f

L(x(t),u(t),t)dt

de zogeheten tobjectfunctionaa1I7. Deze integraal heeft, gezien de rneerdere besturingsrnogeljkheden, ook meerdere

waarden , waarbj uiterard 1(u) L0(t - t0)

Het dynamische optimaliseringsprobleem kan nu gesteld worden als

Breng x(t) vanuit x(tt) ₌ en via

(t) = f(x(t),u(t),t)

door keuze van u zodanig na3r x(t) = Xf dat daarbj de kotn-integraal

t1

1(u)

=

f

L(x(t),u(t),t)dt minimaal wordt."

t o

Voor theorie en toepassingen van Dynamisch Optimaliseren zj verwezen naar de volgende colleges:

a

₇₅

A : Optimaal Besturen van Dynamische Systemen.

(Dr. Ir.J.J.,Kalker)

a 75 B : Nurnerieke methoden voor het oplossen van Optimale Besturingsproblemen. (Dr.C. de Wit) a

₁₆₅

B: Optimaal Stochastisch Navigeren. ( - - )

1.2.

Voorbeelden van Statische Optimaliseringsproblemen.

1.2.1.Constructie van een aluminium luikdeksel met minimaal

gewicht. (Lit.1)

-Q-4. _I

b en t3 zjn gegeven, h en t1 zin "vrij" te kiezen.

1s de keuze van aluminium als materiaal vast ugt,

dan is het gewicht van zo'n luik recht evenredig met

F(t.,h) = h + 120 t

iestric ties

De vervormingsspanning 1, gelijk aan 1800/h,

mag niet

grater zijn dan

_{Ç11 = +50. Dit betekent dat h.}

+ cm.

De buigspanning

= 1500/(t1h) mag een bepaalde waarde

all

niet overschrijden. Dit betekent dat t1h,

6.-f236.

Het buigmoment moet kleiner bl5jven dan het knikmoment.

Dit betekert zoveel als

=O

_{700 t}

_{= Ç}

_d.w.z.

th ? 6.k286

(1+) Tenslotte is er een maximaal toelaatbare doorbuiging,

gelijk aan 1.5 cm. De doorbuiging bedraagt

1+81.71k28

2

t fh

Dit betekent dus dat

_{t1h2> 32l.1285}

(5) Triviaal voor de technicus, maar

van belang voor de,

niet zell denkende, rekenautomaat is voorts het feit

tf. O

Het probleern is dus:

Minirnaliseer F(tf,h) = 120 t.

+ h

onder de voorwaarden

) o

(1)

h -

O

(2)

t1h - 6.Lf286) o

(3)

th - 6»+286

O

()

tfh

_{- 321.l285}

(5)

t1

o

1.2.2. Maximale stationaire stijgsnelheid van een vliegtuig.(Lit.2)

V = sneiheid

V

elevatie van V

= aanvalshoek

rn

massa v/h viiegtuig

g = versneliing v/d

zwaartekracht

= hoek tussen de as

v/d straalmotoren

T

en de nul-liftas

Stationaire toestand betekent krachtenevenwicht

(Horizontaal)

f1(V,,) = T cos(' +

- D - mg sind, = o

(Verticaal)

f2(V,1,û) = T sino

+

+ L - mg cos-= O.

T, L en D zijn gegeven als functies van V en

T = T(V) = stuwkracht van de motor,

L = L(V,.) = liftkracht,

D = D(V,ø) = weerstand.

De stijgsnelheid is

V sin_

De variabelen zijn V,- en

.

Probleem: Maximaliseer V sin

onder de beperkingen

= O

,

f2(YQ,) = O

T

- T(V) > O

1.3. Enige wiskundige notaties, definities en conventies. In dit yak wordt nogal wat vector- en rnatrixrekening bedreven. DaarbU wordt een vector x met regle kentallen

x1 t/rn x opgevat als een (n,1)-rriatrix

Xl

= _{= (x} . . .

. x)

s

X n

X : kolornvector, X r&jvector.

Aangezien x1 t/m x regle getallen zjn, is x

de regle n-dirnensionale vectorruimte.

Vectoren worden vaak zeif van indices voorzien

= (xlk. s s xflk).

BiJ matrices komt de nj-index op de eerste en de kolom--index op de tweede plaats:

fa . . . a

i .11 .ln

A

=1:

(a . . .

. a)

(n,m)

1

nl nfl

De vectorfunctie f heet een re1e m-dirnen5ionale afbeelding van R, als er in rekenvoorsc'nriften

1 t/m

f bestaan, waardoor de rn-vector

in

=

1(&

als beeld van x= (x1

. . . .

x)'

wordt 'oepaald.

Zij f een afbeelding van R in R

Onder de T'afgeleide van f naar x'in een punt _xk,

notatie

_f(xk)

, verstaat men de (n,m)-matnix

(f.

J _{(x )}} i = 1(l)n, j = 1(l)m. L

i

k

element 1e rU en e koloin

Als f een scalaire afbeelding op R' is, dan heet de n-vector

f f

= . s

s

de gradint vanfinhetpunt (.)

. Andere

(b)

iiiI

₌

lxii +

Zij f een scalaire afbeelding op '. Onder de H e s s e

--matrix van f in een punt x verstaat men nu de matrix

Hf(x) ₌

-Ix(&

(12 1)

e]eent

in

de

jer jekolom

Oprnerking: j 2 continu differentieerbare functies is de Hesse-matrix syanetrisch.

Bij het opsporen van een locaal minimumpunt van f in de buurt van een punt wordt vaak uitgegaan van een 2e

orde-ont':;ikkeling van f orn

X0:

f(x) = f(x) + f'(x)(x - + -

)'H()(x

- ) (12.2)

n

= f(0) +>_(x)(xi - x.)

+

- x.)()(x. - x.)

(172b)

i=1

Hierbij is een punt, ergens ophet lijnstuk van x

naar x : =

+ (iÀ)

met O, ) i

De norm van een vector x wordt genoteerd als

Zo'n norm moet aan de volgende eigenschappen voidoen:

llII

\ o

IIXH

= O dan en slechts dan, als x = O

+

z'J<

x

+ ìiII '

ll° II = jol_j

I,

waarbij een re1e scalar is.

Vaak gehanteerde normen zjn

2

(a)

_{UxjI =}

(x1 + . . . .

+ x)

(c) _IlxU = Max (i)

Het punt heet "locaal mininiurnpunt,van f" als

er een ree1 positief getal r bestaat, zo dat uit

ij -

zIL

r voigt dat f(x)

f()

heet "globat1 minimumpunt van f in V" als voor elke

2. _{.inimalisering van functies zonder beperkingen. (Lit.3)}

2.1. Nodige en voldoende voorwaarden voor een vrj locaal minimumpunt.

We beschouwen eerst een _{én-dimensionaa1 geval. Voor}

(0,00) zij f(x) = x + (ln(x))2 . _(13.1)

)I=y, +Y

/

/

/

/

/

Figuur 2.a. Nu voigt: f(x) = 1 + 2 ln(x)/x = 2(1 - _ln(x))/x2

Het nulpunt van f voigt uit _{wton-2aphson-iteratif}

(k+1) (k) x

X

=X

-Voor dit nulpunt x0 vindt men 0.703'+7 (Ga _{dit even nafl}

Voorts blijkt f(x)

_{> O voor ln(x) <1 , d.w.z. voor}

x <e = 2.21828.. Eijgevolg is ±'(x) convex voor x

= X0

Voor elke 4x

met 14xj<2 is

f(x

_{+4x) = f(x) +x f(x)}

_{+ j-x2f (x}

+9x)

= f(x) +

0 _{+ lets,> O} d.w.z. _{is een minimumpunt van f.}

f(x)

= O

Een voldoende voorwaarde is, dat er een getal ¿>0

bestaat, zodanig dat f(x)

_>

O voor elke x

/

x met

Ix -

I

<

, d.w.z. er is een omgeving van x0, waarin de 2e afgeleide beslïst positief is.

In

analogie hiermee kan men uit de formule

t t

f(xx)

=

f(x)

+ x + -jx Hf(

+eAx)

x

aantonen dat een minimumpunt van f is, als voldaan is aan de volgende voorwaarden:

f (x ) = O

x-o

-2e)

Er bestaat een positief getal r, zodanig dat

i,00r alle X x met z - z 'r de

-

-o

-Hesse-matrix Hf(X) definiet positief

is.

Met dit "definiet positie wordt bedoeld, dat voor elke

vector x

/

O het scalaire product

AXHf(X)

x >0 is.

Voor het positief zijn van een symmetrische n3cn-matrix

A =

(a..}

zjn diverse criteria, zoals

D.

_>

O voor i = l(l)n. Hierbij is D1 ₌ Detahk} met h, k = l(l)i

De (reale) eigenwaarden van de (symmetrische) matrix A zjn alle positief.

De aangewezen weg voor het opsporen van een globaal minimumpunt van een functie f ljkt nu:

Zoek alle punten _Xk , waarin f stationair is, d.u.z.

waar

_{= 2}

Bereke voor elk van die punten de Hesse-matrix en

ga na, of deze al dan niet definiet positief is, Eereken in de locale minima de functiewaarcìen. Het punt met de kleinste functiewaarde is het globale

-Helaas is het in de practjk vaak heel rnoeilijk en vooral zeer tjdrovend, orn de elementen van die Hesse-matrix uit te rekenen.

Oak het analytisch berekenen van de gradintvectoren

blijkt in de practijk steeds vaker OD moeiLjkheden te stuiten. Dit kamt, omdt men de diverse minirnaliserings--technieken steeds meer gaat toepassen op functies, warbj het rekenvoorschrift weliswaar énduidig is -anders mag men het geen functie noemen - maar waarbij dit rekenvoorschrift niet is weer te geven door de &n

of andere gesloten analytische vorm.

De in deze paragraaf gegeven beschouwingen zjn dan ook niet voor directe practische toepassingen bedoeld. W&1 is het de bedoeling, hiermee een inzicht te geven in de structuur van 'uncties en in de, op die structuur gebaseerde werking van minimaliseringstechnieken. 2.l.A. Voorbeeld van een

functie

van één variabele, waarvoor

geen expliciete analytische vorm bestaat.

Een gedwongen trillinS kan in cerate instantie worden

beschreveri door de differentiaalvergelijking

y(t) + 5 (t) + 6 y(t) = sin

t

. (15.1) Hierbij is gegeven: y(0) = O en y(7T/2) = i

De algemene oplossing van (15.1) luidt

y(t) = C1e_2t C2e_3t + O.i(sin t - cos t) C1 en C2 volgen uit de randvoorwaarden

0=Ci

_C2

- 0.1 -Tr .-31T/2

- C1e C2e + 0.1

'en vindt C1 = 26.266 en C2

= -

2E.166

Bu nadere beachouwing blijkt de differentia1verge1ijking

jets gecompliceerder te zijn:

y(t)

_{+ 5}

y(t) - y3(t)

+ 6

y(t) = sin t kiiervoor is geen algernene oplossing in expliciete

vorm te geven. Numerieke opiossing met de gegeven

raxidwaarden kan als voigt.

Voer als nieuwe variabele in x Nu ontstaat het simultane systeem

= - 5x

+x3

_-

Gy +sin t

y= X

Hiervan weten we y(0) O en y(/2) = 1 Numeriek oplossen kan alleen als x(0) bekend is. Std x(0) = a Dan wordt x(t) = X(t;a) en y(t) = Y(t;a)

0m te zoren, dat yCIT/2) = i wordt, moeten we a zo kiezen

dat _f(a)

_(y(T172,a)

)2 minimaal wordt.

Deze functie f bestaat en is zeifs 2 continu

diferen--tierbaar naar a. Men kan voor berekening van

f(a)

een numerieke procedure ontwerpen, Met behuip van een

minina1iserings--algoritme, dat alleen gebruik ma&kt van functiewaarden, kan f(a) geminimaliseerd worden.

17

-2.2.

Enige Dractische voorbereidingsmaatregelen.

2.2.1. Schaling van de variabelen.

De grootheden, die bj een optimaliseringsprobleem als

variabelen betrokken zjn, zn vaak geheel verschillend

van aard en dus ook van grootte-orde. Een bruto

scheeps--tonnage van ongeveer 30 000 ton maakt een fout van 30

o

ton acceptabel en bi.j een temperatuur van ongeveer 500

zal niemand schrikken van een fout van 10

Teneinde in een programma voor alle variabelen met een

vrjwe1 gelijke grootte te kunnen werken, verdient het

danrom aanbeveling, van te voren na te gaan, hoe groot

die variabelen ongeveer kunnen zjn.

Heeft een grootheid x bij voorbeeld een verwachtingswaarde

2

c 10' in y en ee-n stap in x van 50 betekent .een stap

in y van 0.01

Ligt de verwachtingswaarde dichtbj nul

,

d.w.z. bestaat

er een redeljke kans, dat Ix/.( <1 zal zjn, dan kan men

beter schalen met behuip van een geschatte bovengrens

x

_sup

van x. Men substitueert dan voor x

x = y

.

(17.2)

Een geljksoortige maatregel kan worden aanbevolen ten

aanzien van partile afgeleiden. Aangezien men bij het

minimaliseren van functies steeds streeft naar zowel

kleinere functiewaarden als ook naar kleinere gradinten,

mag men gevoeglijk de absolute waarden van de

gradint--componenten in het startrunt

als bovengrens beschouwen.

Men bepaalt nu in bet startpunt

eerst het getal

= Max

()j

_.

_(17.3)

(i)

i

De functie f wordt dan vervangen door f

volgens

=

(&/g

= 5000 en wenst men x te bepalen met een tolerantie

van 10, dan kan men de variabele

y

invoeren volgens

y =

_(17.1)

d.w.z. in de functieprocedure wordt

X = y

gesubsti--tueerd. De variabele y varieert dan rond de waarde 1.

De tolerantie van 10 in x wordt dan een tolerantie van

Men kan er in zo'n geval vrij zeker van zijn, dat de

gra--dintcornponenten van f binnen de grenzen -1 en +1

zullen variren.

¿.i.2. Nurneriek differentiren.

In vele optinaliseringsalgoritmen wordt gebruik gemaakt van partile afgeleiden. go'n algoritme bestaat dan in principe uit twee basisprocedures, namelijk

het mininaliseren van f langs een ljn x = x

A2-

en

het selecteren van een zoekrichting in een punt

X.1 ,

dat als rninirnurnpunt op zo'n lijn is geaccepteerd. Wanneer nu zowel bj de 1jnrninirna1isatie als bu het

be--palen van een zoekrichting van gradinten gebruik wordt gemaakt, dan is het bij meer dan drie variabelen in de regel te duur (betr. rekent5jd) orn deze gradinten steeds

weer d.rn.v. numerieke dií'ferentiatie te bepalen.

Is echter de lu,jnminimalisatie zodanig, dat daarbuj alleen

van 'unctiewaarden gebruik wordt gemaakt, dan verdient

het b niet-analytisch gegeven functies aanbeveling,

voor het bepalen van de zoekrichting p in een

start--punt x. de partile afgeleiden te geven d.r.v. een

nurnerieke procedure. i)eze verloopt in principe als voigt.

i1en kiest - over die keuze zell zullen we het hierna nog hebben - een toename-vector

dx = (dx1 . . .

dx )'

-

f

11

De partiele afgeleide

V(x)

an nu worden oenaderd

door n van de twee vo1gede uitdrukkingen

--()

(f(x +

dx.

S.)

-

f(x))/dx. , warbij

(18.1)

L.

=

(&.

. . . .

&.)',

het z.g. voorwaartse

differentiequotint,

!_()

(f(x + dx.

L.)

-

f(x - dx. £.))/(2 dx.) , (18.2) het z.g. centrale differentiequotint

ilet zal duidelijk zijn, dat men in het eerste geval per

gradintberekening n+l functiewaardeberekeningen nodig heeft. 3j centraal differentiren zin dat er

19

-Daar teenover staat, dat de fout in de

benadering van

bij gebruik van centrale differenties

_{voor naar nul}

da±ende dx

_{sneller naar nul convergeert dan bij}

_voorwaarts

nur.ieriek differentiren.

Aangezien evenwel de partile afgeleiden alleen

maar

hulpmiddelen zijn, die kunnen dienen tot het snel

_opsporen

van het minirnumpunt van een functie, telt dit

nauwkeu--righeidsargurnent minder zwar dan het feit, dat

men

bij centraal differentiren ruwweg twee keer zoveel

functiewaardeberekenìngen nodig heeft als bj gebruik

van voorwaartse differenties.

Voorwaarde is echter wel, dat men bij numerieke

gradint--berekeningen de al dan niet voldoende benadering

van

het niinimumpunt nooit mag laten afhangen

_{van het al of}

niet voldoende klein zi.jn van de

- benaderde en dus

mo--ge]Jjk foute - gradintcomponenten.

en aanzien van de keuze van dx

_{kan het volgende}

wor--den opgernerkt. i3eschouw een functie f

van

n

varia--bele x. f zj minstens twee rnaal continu

_{differentieerbaar.}

Dan geldt:

(f(x + h)

_{- f(x))/h = f'(x) + 3 h f(x + eh)}

_.

_(19.1)

De afbreekfout in deze benadering is dus

Daarnaast is er een computerfout in f(x+h)

en

die in de benadering van f(x)

_{een fout kunnen}

veroor--zaken, die in absolute zin wordt gemajoreerd door

_{l0°jf(x)I/h,}

;.angezien bu afnemende h de afbreekfout daalt

en de

computerfout stijgt, is de beste waarde

van h de waarde

h°, waarbij deze

_{twee fouten even groot zijn. Deze h°}

voigt dus uit

3h0jfu(.)

=

h°

= l05(21f(.)/f"(.)I)

_(19.2)

Voor f"(.) kan men b

benadering nemen

(f(xh)

- 2 f(x) +

f(X_h))/ha

_{met een rede-}

_(19.3)

Deze bealing van h° kan, als de

2e

afgeleide wat klein

uitvalt, leiden tot een onwaarschijnlijk grote waarde van

o

h

.

angezien men ecnter niet zeker is, dat deze

ai--geleide oak elders zo klein is, kan men

,

als de met

(19.3) en (19.2) berekende h° 2-roter dan 0.001 blijkt

te zijn, voor h° = 0.001 nemen.

Voor de bepaling van dx kan men dus werken met de

vol--gende regels

Voor i

:= l(l)n

Bepaal een benadering van -(.) uit

xi

= (f(

+

/1000) + f(x

-

g/lOOO)

-

2*f(x))1O6

i

dx.

10

(2*f(.)/ --(.)D

ax.

als dx.> 0.001

,

dan

dx

0.001

Zie ook Lit. 4.

(20.1)

(20.2)

21

-2.3. Ljnminima1isaties. (Lit. ₃₎

Bij diverse rriînimaliseringsmethoden is het minimaliseren

van de gegeven functie f langs een bepa1de Ljn een steeds terugkerende subroutine. Uitgaand van een

start--punt wordt een zoekrichting bepaald, waarna f

langs de L,jn x

= X1 + .

_{wordt geminirrialiseerd. Op}

die Lin wordt f(x) een functie van X

f(x

+X.) =

f(X)

(21.1)

Bij een z.g. kwadratische Ljnzoeker gebeurt dit minima--liseren uitsluitend met behuip van functiewaarden. iIen onderscheìdt twee fasen, n.1. de stapzoekfase en de locale minimalisatie. In de stapzoekfase wordt

stapsgewi,js veranderd, d.w.z. ,> O, s, 2s, 3s,

totdat men in de situatie is gekomen, dat voor de

punten a ks met

_f(a)

= fa b = (k+1)s met f(b) = fb en

c =

(k2)s met f(c)

= fc

geldt, dat zowel fb fa als fb fc

Hoewel de zoekrichting steeds zo wordt bepaald, dat f in die richting vanuit dalend is, kah bet voorkomen, dat reeds voor a = O en b = s de functiewaarde ft

gro--ter is dan fa. In dat geval wordt s telkens gehalveerd, totdat ook dan een "dalsituatie' is gevonden. Aan het eind van de stapzoekfase ugt in elk geval het minimum--punt van f op het interval [a,c]

Men legt nu door de punten (a,fa), (b,fb) en (c,fc)

eeri dalparabool, d.w.z. een kwadratische functie

qu(x)

= px2 + + r voldoet aan

pa qa+r=fa

pb2 + qb + r = fb

PC2 + qc 4 r fc

Van deze functie is het minimumpunt

.

/

-, a (ib - fc) + h (Ic - fa) + c (fa - fb)

- q,p

(fb - fc) + b(fc - fa) + c(fa - fb)

Voor dit punt d wordt nu fd = f(d) berekend.

Indien b1jkt, dat d voldoende dicht bi.j b ugt, dan wordt d als minimumpunt van f _{geaccepteerd. I} dat niet za,

Indien berekening van de afgeleide van fk naar À weinig moeite kost, kan men dez.e informatie evenzeer rebruiken. Men kan dan met de stapfase doorgaan, totdat men twee

punten a en b heeft gevonden, waarin dfa = fr(a) < O en dÍ'b = f(b) > O is. Met de gegevens a,fa,df:,b,fb,

dÍ'b kan men f (À) op ca, benaaeren met een . graads

polynoom. Hiervan is het minimumpunt a op eenvoudige wjze te berekenen.

Bljkt dfd = f(d) voldoende klein, dan wordt d als minimumpunt van f geaccepteerd. Is dat niet het geval, dn kan men a of b door d vervangen, al naar gelang d±'d negatief dan wel positief is.

Betreffende de keuze van s verschillen de gebruiken nogal. Wenst men een staplengte van rn

cordinaateen--heder. - men werkt meestal met m = i , hoewel bj geschaalde

cordinaten m.= 0.1 logischer lijkt - dan moet men voor

s nemen: h = m/I),Ij met = . (22.1)

j=l

Vaak kent men de richtingsafgeleide van f langs de ln in het startpunt X1

f(o)

=

Met een schatting f van het minimum van f en een

kwa-est

-dratisci-i verloop van f(X) nabij het minimumpunt is de

stap naar het minimum in cordinaateenheden:

mt = 2(f - f(x.))/f(0)

est

i

De À-stap zou dan worden:

k = 2(f - f(x ))/pf (x )

est

o

i X O

Voor s neemt men dan de kleinste van de wa7rden h en k.

23

-2.1+. _{Minimaliseringsmethoden, waarbij geen gebruik wordt}

gemaakt van gradinten.

2,--.l. De Directe Zoekmethode van Hooke en Jeeves.(Ljt,

₅₎

Deze methode bestaat uit een locale zoekprocedure E(s1,x) - exploratory move - en een gerichte bewegingsprocedure. In de zoekotapprocedure wordt vanuit een punt x docr successieve cordinaalwijzigingen van + gezocht naar een Dunt met een kleinere functiewaarde. Zie verder het stroomdiagram op biz. 2+

De gerichte beweging zie stroomdiagram op biz. 25 -werkt zo, dat na een succesvolle zoekstap de beweging in die richting wordt voortgezet en bj herhaald suoces

versneld.

B1,jkt het lopende punt x echter na zo'n gerichte bewe--gins een grotere £unctiewaarde op te leyeren dan het, na de laatste zoekbeweging bereikte punt y, dan wordt

X weer naar teruggezet.

Heeft een zoekstap geen resultaat, dan wordt LX verkleind.

Heeft een zoekstap mt

<geen succes, dan wordt

het bereikte punt als minimumpunt geaccepteerd. De fout in de cordinaten van dit "minimumpunt" is dus , als

de methode goed werkt, kleiner dan

L.

De gerichte bewegingen en zoekstappen worden uitsluitend door kwalitatieve veranderingen van f(.) bepaald. De cordinatenstap is voor alle cordinaten even groot. Ilet is dus dringend gewenst, de variabelen te schalen, zoals aangegeven in 2.2.1.

Na bet stroomdiaram op biz. 25 staat een uitgewerkt

voorbeeld. i{ierb5j is

f(x) = 52

4

_{+ 72 x1x2 + 73}

_{- »+8}

_{- 36+}

_{x2 (23.1)} Startpunt (5,3) . Aan het begin is = 0.5 . Later

wordt telkens met een factor 0.2 verkleind. Na een succesloze zoekstapprocedure met een _{<O.0l breekt} de algoritme af.

inarray(x) si k xn X fx si Xflk

Xflk + C

fxn Ck := - Ck Xfl1 _Xflk + 2c1 fxn

f()

Xflk

Xflk -

Ck k k L3 fx: frs 'y.

k

si fr end fx fxn := _Xflk

'I,

gin

input:

r;

x; sO:=sl:=f(x)

'±or' i:=i()n 'do'

c.

¿:= r

'for'

i:=1(l)n

'do' c.

_:=L\

-

25

-Directe zoekmethoäe

end

z := := X

:= 2x - z

_; sO

:= si

si

:

nee

si := sO

.12 _{SO :=s}

2

(k.5,2.5)

29.25

5 3

(+,2)

- 20

(k.5,2.5)

293,;:5

(5,3)

6 k

-20

5

(.5,i.5) -23k.75

8

(2.s,O.5) -363.75 (3.5,1.5) -23k.75 (k.5,2.5)

9 7

-368.75

8

(2,1)-1+35

12 9

(0.5,0.5)-256.75

(2,1)

-1+35

(3.5,1.5)

13

10

(2,1)

_L35 11

(2,1.5)

k53.?5

16 12

(2,2)

-2+T6

(2,1.5)

-1+53.75

_(2,1)

17 13

(2,1.5)

-1+53.75

1k

_(1.5,2)

Lf75 ₂₀ 15

(1,2.5)

-1+69.75

(1.5,2)

k75

_(2,1.5)

21 16

(1.5,2)

-1+75 17

(1,2)

-1+88 2k 18

(0.5,2)

-1+75

(1,2)

-1+38

(1.5,)

25

19

(1,2)

-1+88

20

(1,2)

-1+88 29 21

0.1

22

(1,4)

-1+88 33 23

0.02

2k

(1,2)

-1+88 37

25

0.001+.001+ 26

(1,2)

2+88 26

(1,2)

2+88 1+11+1

27

0.00k<0.01 ,

afgelopen, uit

27

0.00k<0.01 ,

afgelopen, uit

27

--.. Neider & Meads Simplex Algoritme. (Lit.

₆₎

Wanneer functiewaarden op de n of andere fysische manier worden geproàuceerd, dan zjn deze waarden door--gaans behept met ruis. Nu zijn er wel filtertechnieken, die die ruis er zo goed moeljk uit kunnen filteren, maar echt deterministische, van alle ruis ontdane,

functiewaarden krjgt men nooit. Als men met zlke functiewaarden werkt, kan men het begrip afgeleide wel vergeten. Zeifs de aanname van een locaal convexe

functie is dan vaak niet eens verantwoord.

Voor het minimaliseren van zon functie is door Neider en Mead een techniek bedqcht, die de situatie telkens in een "clusters' van nl punten beldjkt. De algoritme werkt als voigt.

Zj x een n-vector en f een scalaire functie van x.

Men kiest een 'simplex" (driehoek in R2, vierhoek in R-p)

van n. punten

x t/m x. Men berekent voor

i = 0(i)n: y1

= f(x.)

. Van deze punten bepaalt men het maximumpunt x, en het rninimumpunt x1 , d.w.z.

voor alle

-h en

1)

L

Cok bepaalt men

/n en =

f(s)

ih

Deze simplex wordt nu gewijzigd, d.u.z. men vervangt in eerste instantie &n punt,

h door een ander punt

met een kleinere functiewnarde. Hiervoor zjn drie standaardbewerkingen met de stoere namen

reflectie d.m.v. een parameter

_>

O

expanie d.m.v. een parameter _{.->i en} contractie d.m.v. een parameter

_fi

(0,1)

Allereerst wordt Xli gereflecteerd naar _x voigens

= (i01)E

y = x

Indienu Y<YL blîjkt, dan wordt dit

verder gxpandeerd;

=

'??

-

, (27.5)

=

f()

-

-B1jkt nu, dat y

ç

, dan wort x. acor x vervangen. In het andere geval wordt Xh door x vervangen.

Is ?yL

, dan wordt nagegaan, of y> y1 voor alle ih.. Is dit niet zo, dan wordt x aisnog aan de simplex toe--gevoegd als opvolger van Xh

Blijkt y>

h' dan wordt contractie toegepast . Bl5jkt

y wel > y. voor i / h , maar

h dan wordt

door x vervangen. Ook nu wordt gecontraheerd;

=flh +

(1

L

ßljkt nu =

>

h dan vindt totale

contrac--tie plaats, d.w.z. voor alle i wordt

Is rd

h vervangen door

Voordat een nieuwe simplex wordt bewerkt, vindt een convergentietest plaats. Deze is statistisch van aard.

?ien gaat n.l. na of de standaardcleviatie voldoet aan

-2

- y) /n 1=0

<cE..

Op de hierna volgende bladzjde staat een stroomdiagrarn van deze minimaliseringsmethode.

(28.1)

(28.2)

reflectie

x:(1 +o) -dx,;

y f(x); nee 'V, X y

:=x

: y tractie :: := +(i -y : f(x ); inarray ;

Minmax(xy,xi,yi,xy);

e gi r X. (x. + x -1

-i

1

y1 f(1); nee expansie := X

n

y geconvergeerd. end.

fase.

2»4.3. Kwadratisch convergente ininimaliseringsmethode van Powell. (bit. 7) Powell's methode rnaakt alleen gebruik van functiewaarden. In tegen--stelling tot }Iooke & Jeeves' en Neider & Mead's methoden wordt hier echter van de functie in de buurt van het minimumpunt een bi benadering kwadratische structuur verwacht.

Bij een zuiver kwadratische structuur en natuurlijk voor een convexe functie met een definiet positieve Hesse-matrix convergeert Powell's

methode in n iteratiestappen naar bet minimunipunt. Daarbij bestaat elke iteratiestap uit n+1 li.jnminimalisaties. Steit men voor elke

lijnminimalisatie het benodigde aantal functiewaardeberekeningen op q dan wordt het kwadratisch minimum dus bereikt na qn(n+1) functiewaardeberekeningen. Dit is vrij veel. Bij kwadratisch con--vergente methoden, die wi van gradinten gebruik maken, bestaat elke iteratiesta uit

een numerieke gradintberekening en

een kwadratische lijriminimalisatie.

Dit betekent, dat zo'n methode per iteratiestap n(n+a) f-waarde--berekeningen vergt.

Reeds bij n=3 en q5 is Powell's methode al in het nadeel ten

oD--zichte van zo'n kwadratische gradintniethode:

Powell's methode : 60 f-waarden per iteratiestap,

Kwadratisch cony. gradintmeth. : 21 - " - " - " - " -.

In feite is er dus weinig reden voor een uitgebreide behandeling

van deze Powell-methode.

Voor een uiteenzetting betreffende bet principe van kwadratische convefgentie etc. wordt verwezen naar 2.5.

Korte beschrijving van Powell's methode:

En iteratiestap heeft k fasen. De ke fase bestaat uit drie delen:

(i) Met een startpunt en zoekrichtingen t/m wordt voor

i = i(i)n bet lijnminimum x. van f beiaald op de lijn

/ \

-k)

x=x.

± D.

-

i--1

(2) Men kiest nu

_{n+1 n o}

= x - x en bepaalt het minimumpunt van f

op de lin x = x k k Dit is dan het startount voor de

(k+l) (k)

men

=

(k+l ) (k+1 De procedure breekt af, ais!1 x

_n

- x

_o

Het eerste startpt (1) _{wordt door de gebruiker opgegeven, Voor}

de zoekrichtingen in de eerste fase neemt men

D)

= d.w.z. p = 31 (3) Van de zoekrichtingen u (k) (k) afstand

fix.

-

x. (k) . wordt het

diegene gescbrapt, waarvoor de

s

rootst was. De daarna volgende t.

-i

-i i

2.5. Gradi nt-ininimaliseringsniethoden.

2.5.1. De anti-gradint of steilste-dalingsmethode.

Bij de anti-gradintmethode moet men beschikken over numerieke of analytische procedures ter berekening van de componenten

van de gradintvector f(x).

f(x) =

richting van de steilste daling

In bet punt x kan men in alle mogelijke richtingen infinitesi--male stappen van lengte ds nemen. Definieert men zo'n richting door middel van de eenheidsvector

a=(a

. . .

.a)

met

a.1

- i n i

1=1

dan wordt de infinitesimale toename van f over zo'n stukje n

df(x ,a) =

o

x.o i

(x )a.ds = af (x

- xo

) ds

i

-i=i

Krachtens de eigenschapren van Euclidische inproducten is nu df(x ,a) als functie van a minimaal, als men voor a een eenheids-

-o--vector neemt, die tegen de gradint in is gericht;

=

-(32.1)

(32.2) Opgave: Gegeven is een functie g van de variabelen a1 t/ra a1

g(..) =f.a. + f(l

- (32.3)

Ees dat g(..) minimaal is voor a.

- f./(f), i1(i)n-1

Wanneer men bet niinimumpunt van f(x) zoekt, ugt het dus enigszins voor de hand, dat men vanuit een stetpunt X gaat zoeken in e

richting van de anti-gradint. Er zn nu in principe twee klassen

van steilste-dalingsmethoden:

(i) Men start in x

_o

en ininimaliseert f(x) langs de lijn x = x -Af (x ).

-

- o

men zoekt van hieruit weer in de richting van de antigradint, nu berekend in x.r Dit herhaalt men, totdat de iteratie af--breekt op grond van een afbreekcriterium, zoals

- < c of

'k

- k-1 <

Deze zo genaamde "lange stap gradintmethode" is in de

prac--tijk weinig succesvol gebleken0

(ii) Vanuit het startpunt construeert men de numerieke oplossings--kromme van het stelsel van differentîaalvergeljjkingen

dx.

i 3f

- - X i-

i = 1(1)n,

De groote van X kan men dan nog variren, afhankeljjk van de norm van de gradint van f0

Oak deze methode heeft vele practische bezwaren. Nabij het

minimumpunt bljjkt de numerieke oplossing ?f vroegtijdig

vast te lopen f eindeloos heen en weer over het minimumpunt

te schieten.

2.5.2. De geconjueerde gradintrnethode. (Lit.

_{3 & 8)}

Deze methode is oorspronkelijk (1952) ontworpen door Hestenes &

Stiefel, is in 1964 als optimaliseringsprocedure geprograinmeerd

door

Fletcher & Reeves

_{, terwiji o.rn. door Kowalik in het zeer recente}

verleden hierop enige belangrijke modificaties zijn aangebracht.

Men is oîj het ontwerpen van deze methode van twee principes uitgegaan.

Allereerst wordt f(x) in de buurt van het minimumpunt 'benaderd door

een

2e

_{orde Taylor-ontwikkeling met een constante positief-definiete}

Hesse-matrix:

f(x) = f(0) + f()x + xH(ex)

-c ± bx

T

xA x =f(x)

In de tweede plaats worden van deze kwadratische 'benadering f

van

f een n-tal onderling "geconjugeerde" richtingen

D

t/m

bepaald.

Voor twee geconjugeerde richtingen

D. en D.

geldt per definitie:

n. A n.

,

d.w.z.

_-1

j

i

ij

£ A

.

= O voor i

j

Men noenit die richtingen ook

Ud

A-loodrecht. Ret grote voordeel

van die richtingen is het volgend.e:

Bij een willekeurig startpunt x

kan men voor i = 1(1)n de lijnminima

x. bepalen van fop de lijnen

X = ₊ .

Als nu

t/m

_£nl

onderling geconjin.geerd zijn, dan is

X

het niinimumpunt van fq

Deze eigenschap zal met een 2-dimensionaal voor'beeld worden toegelicht.

Zij f (x) = x2 + x x

+ X2

-

9x

- 9x

1

12

2 1 2

=

,(

+ (-9

-9)

De niveaulijnen van fq() zijn ellipsen met x

= (3

3)

als

rniddel--punt.

(3I.1)

- 3

-We kiezen als eerste zoekrichting = (i o) en als startpunt

x = (o o) . Wanneer we nu f(x) minimaliseren langs de lijn

x=x.

-

-o -u

₀₂

dan wordt op die lijn f(x) = f(À) = - 9''.

f*

is minimaal voor = 14.5 en het eerste lijnminimum wordt

=(14.5 0)

Dit punt is het midden van de koorde, die de niveaulijn f(x) = O van de X1-as afsnijdt. Het is ook het raakpunt van deze X1-as met de niveau-ellips door dat punt. De, aan de X1-as toegevoegde richting

is de richting van de middellijn door het punt (14.5 o). Voor die

richting geldt :.

,(2 _i'(i\

= O . Hieraan voldoet _{(-i 2)}

1

2/0)

Minimalisering van f langs de lijn x = X1

+ , d.w.z.

,zodat

f(x) = 3) - 9..\ - 20.25 , levert als minimumpunt

= 1.5

' 2 =

Bij de geconugeerde gradintmethode worden dq geconugeerde rich--tingen bepaald zonder evaluatie van de locale Hesse-matrix.

Heeft men n.l. het minimumpunt x. van f op de lijn x=X. + À1.

1

-

-i--1 -i-1

gevonden, dan kan een, aan _i1 geconjugeerde richting worden gevonden door de gradinten g. = f (x. ) en g. = (x.) te

be--j.-1 i-i-1 i

-rekenen. Voor D. kan men dan nenien 1

+ 7

£i_1.i_1

Voor de eerste zoekrichting neemt men meestal p-g=- f (x

-o -o

x-o.

Het complete rekenschema voor het verrichten van n iteratiestap

bij bet mininialiseren van een willekeurige, in het algemeen niet

kwadratische functie f van x luidt nu:

(i) Beginnend bij een startpunt x berekent men = f (x ) en D

-Voor i = 1(1)n bepaalt men bet minimumpunt x. van f op de lijn

x=x.

+>.

-

-i-1 -i-1

In bet gevonden minimumpunt X. berekent men

= en daarna de nieuwe zoekrichting = - + met = D--i-1 _(35.1)

Blijkt na een lijnminimalisatie, dat in het bereikte

punt

x. de

gra--dint voldoende klein is, cLw.z.

k±k

dan stopt de algoritme. Het

punt

x. wordt dan als minimumpunt van

f geaccepteerd.

Blikt na n iteratie, d.w.z. na n lijnminimalisaties, het bereikte

punt x volgens dit criterium nag niet acceptabel, dan neemt men x als nieuw startpunt x . De cyclus van lijnminimalisaties en

be--n

ü

-paling van zoekrichtingen begint dan overnieuw.

Als convergentiecriterium wordt ook vaak de uitkomst van de vraag 1((k) (k)Jf

<

?

(k) (k)

gehanteerd. Hierbi zijn x en X respectieveli.jk het

begin-en eindpunt van de ke iteratiestap. Men heeft n.l. vaak geconstateerd, dat er een onzuiverheid in de gradintberekening kan optreden,

waardoor de gradintcomponenten nag niet voldoende klein lijken,

(k)

terwijl

hec lopende punt x

practisch niet meer verandert.

Dit

kain o.In. vooi- als de

gradiëntberekenirig gebeurt d.m.v. numerieke

differentiatie. Het verdient dan oak aanbeveling, de berekening van de toenamevector dx (zie biz. 20) bij het bereiken van het nieuwe

(k1) (k)

startpunt x =x nag eens te herhalen indien b.v. blijkt,

-a

-n

dat

ç(k) _x(k)ff>o.5

Voor niet-kwadratische functies bestaat de mogelijkheid, dat de

aan--naine " f(x) f(x) gewoon nergens op slaat. Dit is o.m.

merk--haar, als men een zoekricl-iting vindt, waarvoor f(x.

+ À)

als

functie van A niet af-

maar toeneenit, wat

niet bepaald de bedoeling kan zijn. Teneinde dergelijke absurditeiten, die zo'n kwadratische methce van een zegen tot een vloek maken, te elimineren, heeft o.in.

Kowalik het volgende bedacht.

37

-De functie 1' moet in de nieuw beDaalde zoekrichting in elk geval

dalen. Dit betekent, dat de projectie van . op

richting van

.

moet hebben. Het inproduct p.. is Bovendien eist men dan, dat die proectie minstens

O.2 j moet hebben, met j =

)2

Dit

dat voldaan moet zijn aan de voorwaarde

- O.2*j p.

j*j

. I

-i

.iII

Blijkt dit niet het geval te zijn, dan wordt

. =

-genoriien.

Hoewel dit z.g. "resetten" nog niet voldoende is getest, heeft Kowalik reeds in enkele gevallen een belangrijke bekorting van de rekentijd

g e constate er d.

de tegengestelde

dan negatief.

een lengte van kamt erop fleer,

2.5.3. Methode van Newton.

Meer ter voorbereiding van het principe van de Quasi-Newton methoden

dan om het belang

van

de methode zelf,volgt hier een korte

behandeling

van de minimaliseringsmethode van Newton.

We beschouwen allereerst een functie f van gn variabele X. Het miniinu.mpunt van f zij gelocaliseerd op een interval [a,b], d.w.z.

we weten zeker dat

er minstens n punt c [a,b] is waarvoor f(c) f(a) en f(c)

f(b).

f"(x) > O op [a,bl.

Het tweede gegeven vertelt ans, dat de functie ap [a,b] convex is,

d.w.z. de eerste afgeleide neenit monotoon toe en de grafiek is 'glad"

met de bolle kant naar beneden.

J I

a b c

Er is dus precies n [a,bl met f'(x) = O.

x is het minimumpunt

van

f op [a,b].

Nu is f'(x)

= f'(a) + (-a)f"(),

a < <

Met f

T()

= O volgt nu, dat

X = a

-Het beroerde is nu, dat we niet kennen.

We itereren nu naar x volgens:

- (o)

(ji)

(k+1) _: -

f'(x)/f"(x)

Voorbeeld: f(x) = X3 + X2 -- 2x

+ 5

f(0) = f(1) _{= 5} f(O.5) = .375

f'(x) =

3x2 + 2x - 2 f"(x) = 6x + 2. = O voor

= (/T-)/3 =

O.5858371O.

(38.

1) (38.2) (38.3)

-39-Iteraten van : (volgens (38.3))

k (k) o

0.585837704

Afschrikwekkend voorbeeld.: f(x) = cos x - 0.01 (k+1) (k) (0)

Opgave: Verklaar de foute resultaten bij de startpunten x =

0.8

en

0.9.

Analoog aan het an-dimensionale geval rnoet men voor X uit kunnen gaan van de volgende vooronderstellingen:

Er is een gesloten en begrensd eied U c waarbinnen de functie f

in minstens n punt kleiner is dan aan de rand.

De Hesse-matrix H(..) is positief definiet voor elke E

J.

In dat geval is er precies n rninimumpunt van f. Nu geldt f (x) = f (x

x

) + H ()(-x ).

x-0

f

0

(39.1)

Iteratie: X = x (k) K (k)

f(X/ftt(X)

0

0.8

0.9

0.98

1 o,691

2,706

0.9935

2

0.150

_3.176

_0.9906

3

0.008

3.1i6

_0.9886

0.010

3.116

_0.98805

5

0.010

0.98802

6

_0.98802

1

0.625

2

_0.5516

3

o.5L86

ln(1-x). f'(x) = - sin x + 0.01/(1-x) f"(x) = - cas x +

0.01/(1-x)2

f'(x) =

O voor x =

_0.98802

Met f(x) = O voigt nu, dat

= -

H1()f(x0).

x kan iteratief worden benaderd. d.m.v.

(k+1) (k)

H1

(k) (k)

(ii) x = x - (x )f (x ). (IO.2)

-

f

(k)

Als goed is, d,w.z. ais Hf(x0) > O is, dan kan x bijzonder snel naar x convergeren. Vaak bestaat evenwel het gevaar, dat de stap

(k) (k+1) (k+i)

van x naar x voert, waarbi X buten het gebied U valt.

Het iteratieproces convergeert dan rnisschien nag wel naar een

minimum-punt, maar zeer zeker niet naar x U.

(X)

I

f

?X

e Neon-methoe heeft nog meer nadelen. Zo moet in elk punt (k) de Hesse-matrix worden berekend en geinverteerd. Deze informatie wordt bÍj het voigende punt in het geheei niet meer gebruikt.2

Vooral wanneer het berekenen van de partiëie 2e afgeleiden cJ J.

i

J

iangs n'merieke weg moet gebeuren is dit een (reken)tidrovende zaak, vaardoor het voordeel van de (mogeiijk) seflere convergentie weralt tegen de nadeien, dat

men weinig of geen contrie heeft op het "wegschieten" naar een ander extreem

het aantal f-berekeningen voor numerieke benadering van de elementen van de Hesse-matrix ruwweg Bx za groat is ais voor numerieke benadering van de gradientvector d.m.v. voorwaartse different ies.

2.5.4.

Quasi-Newton methoden.

Wanneer men een kwadratische functie

f(x) = x'Ax + b'x + c

minimaliseert d.m.v. lijnminimaiisaties van f op de lijnen

x=x.+X, i=0(1)n-1,

dan kan de berekening van de gradinten g0, g1 etc. een huipmiddel zijn tot het iteratief benaderen van de inverse van de Hesse-matrix A

van f d,na.v. matrices H0, H1, H2 enz.

Hiermee kunnen dan de zoekrichtingen p0, etc. worden bepaaid volgens

ID. = - H.

i

Deze richtingen blijken dan onderling A-ioodrecht te zijn,

d,w.z. p! A p. = A. 6.

-i

J

i 3J

De voornaaznste q.uasi-Newton methoden zijn die van Davidon, Fletcher & Powell en die van Broyden. In de zachte-waren-bilDiiotheek van het

T.H.-Rekencentru zijn vier voorgerrogrammeerde procedures beschikbaar, d0w.z. aanroepbaar in Algol en Fortran. Dit zijn

drie D.Fl,P-versies:

optflp, optfpo, optfpn en een Broydn_versie: optbro.

De Davidon-Fletcher-Powell methode werkt in principe als voigt:

Men kiest een startpunt

De matrix H0 wordt als gnheidsmatrix genomen.

Men berekent g0 =

Men bepaalt voor i = 0(1)n-1 het minimumpunt x.1 van f op de iijn

X = X. + X0.

-

-J_ z_i

Hierbi,,j is p. = - H.g.

-.1 i-3.

3. x.1 is uitgangspunt voor de volgende linzoeker.

De matrix H.+1 wordt ais voigt geconstrueerd:

s. =x.

-X.

-i

i+1

-i

= i+1

-(IL1.2)

A. =

B.

-(H.)(H.yj'

H.

= H. + A. + B..

j.+1

i

J.

i

p .yH

H. y. p

_iJ-i

-t t

ii

ii

5.

i:i+1. Ga terug naar 2.

(L3 .6)

Als het niinimumpunt na n stappen niet is bereikt, dan gaat raen

gewoon door. Soms wordt voor i

n, 2n, 3n enz. volgens bovenvermelde

formules verder gerekend, soras ook wordt de H-matrix weer

terug-gezet naar de eenheidsmatrix.

Dit gebeurt trouwens ook als blijkt dat

_{i__i}

!g.

> O is, zodat de

_.-.

zoekrichting p. niet in de richting van kleinere functiewaarden wijst.

Men zij er dus voor gewaarschuwd, dat de H-matrix bij het bereiken

van het rainimumnunt in het algemeen geen goede benadering is voor de

inverse Hessiaan van f.

Het miniraumount wordt "bereiktt' als voldaan is aan het

n of andere

convergentiecriterium.

Broyden's methode wijkt alleen af van die van Davidon c.s. voor wat

betreft de iteratiestappen in de berekeningen van H1, H2 etc.

De methode werkt als voigt:

Schat

en neem H0 = I(n,n).

Bereken g

_-o

= f (x ).

i:0.

Bereken

i

p. = - H.g.

_iJ.

en minirnaliseer (ten aanzien van x) de functie

f(X) = f(x. +

Dit minimum wordt bereikt voor X = À..

_i

Neeinnux.

_-2+1

₁

x- +X.n..

_it-i

Bereken

_{£i+1 =}

en ga na, of

g.1

_{< c.}

Is dit zo, dan stot de algoritme. Is aan dit convergentiecriterium

niet voldaan, bereken dan:

= .i+1

-1ii

-i

H.

H.+

-r

2.5.5.

Kort

overzicht

van

bet, op bet T.H.-Rekencentrum beschik'bare,softwarepakket.

10.301: optdzm (x, n, fx,csteD ,r , delta, f, cony,

meval);

Dit is een Hooke & Jeeves-versie, waarin, na een 'pattern move' van naar x, niet wordt nagegaan of f(x)

f(y).

Hierdoor ontstaan overbodige 'exploratory searches'.

Bu een delta,die ugt tussen 0.01*cstep en 0.1*cstep, geeft

deze methode een snelle en doorgaans goede startwaarde voor een TTfijnere'

rninimaliseringsinethode.

Waarschuwingen:

w1: schaal vooral alle variabelen op waarden rond 1.

w2. Het gevonden "minimumpunt" is beslist onbetrouwbaar, als de gradiënten dicht erbij erg groot kunnen zijn.

10.301: optpow

10.308: otplo

Twee versies van Powell's gradiëntloze methode, beschreven in 2.1.3. optpow bevat een aantal programmeerfouten, waardoor de procedure in volkomen reguliere gevallen stuk loopt.

optplo is eenrersie van Loots±a, die het, o.m. blijkens de ervaringen van de Centrale WerkgroepWiskunde van de Afd.- Scheesbouwkunde,

uitstekend doet. De Kwadratische lijnzoekers breken snel af, zodat het aantal f-berekeningen niet al te groot wordt. Doordat gradinten niet worden berekend, kan er ook geen fout optreden als gevolg van een

foute nu.merieke benadering van die gradiënten.

10.302: optflp

10.303: optfpo

Twee versies

van

de quasi-Netonmethode van Davidon, Fletcher & Powell. Bij beide versies worden in elk punt x zowel de functiewaarde als de gradì nt berekend.

Optflp werkt met een kubische li,nzoeker, optfpo met een kwadratische lijnzoeker. Optflp werkt dan ook jets sneller.

Beide methoden zijn alleen toepasbaar op functies, waarvan de gradient analytisch leverbaar is.

Afbreekcriteria:

Optflp breekt af, als

ann de volgende beide eisen is voldaan:

(i) Voor i > n is de lijnstap vanuit het startpunt X. tot aan het lijnminimum

x. kleiner dan c:

1+i

lx.

-x.

H

<E.

1

i+1

'Eucl.

(2) <

OptÍpo breekt af op hetzelfde criterium. Alleen wordt hier als norm van "r

xehanteerd: ¡jxj

10.304: Optcgr

Dit is een versie van de econ,ugeerde gradintmethode van Fletcher & Reeves. Ûok hier worden in elk punt x zowel f(x) als g(x) = f(x) 'cerekend.

Oak deze methode is alleen 'oruikbaar, als de gradient analytisch

lever-baar is. De methode convergeert duidelijk minder snel dan optflp, hoewel ook hier de lijnzoeker met kubische interpolatie werkt.

N ₂

De procedure breekt af, als g' = E g. < _E.

i 1

Wanneer men dus "oi,j otflp met c = iO werkt, moet men hier voor eenzelfde nauwkeurigheìd met c = 10 werken.

10.306: Optbro.

Dit is een versie van Broyden's quasi-Nebonmethode.

Bij deze methode wc:den de functiewaarden en de gradiënten door twee aparte procedures geleverd. De linzoeker is kwath'atisch en werkt dus alleen met functiewaarden. Hierdoor is het mogelijk, de gradiënten

ook via een numerieke procedure te berekenen. Alleen nioet men dan

voor "eens en altijd" de toenamevector dx vaststellen. Ook deze methode convergeert,als

g'g < .

Berekent men de gradiënten numeriek, dan is dit criteriurn mogelik minder bet rouwbaar.

10.305: Optfpn.

In deze Fletcher & Powell-versie wordt de gradiëritvector aan het begin

van elke lijnminimalisatie numeriek bepaald. In het startpunt

worden voorwaartse differentiequotiënten bepaald met behulp van een, door de gebruiker opgegeven vector dx, d.w.z. voor i = i(i)n neemt

men = (f(

+ dx.)

-Na elke lijnminimalisatie wordt deze vector dx dan za goed mogelijk gekozen, gebruik makend van Stewart's methode. De lijnzoeker werkt met kwadratische interrolatie, waarbij alleen ñinctiewaarden en de

145

-De procedure rekent o een functie met nul als minimumwaarde. Wanneer men dus een schatting f . van dat minimum heeft en f . is negatief,

min min

dan moet men f(x) vervangen door f(x) - f.. Men moet er daarbij zeker van zijn, dat f(x) > O blijft.

Wil men b.v. f(x)

=

h.(x)

minimaliseren (zander restricties t.a.v. x= (x1...x)') dan vervange men f(x) door f(x) +

Het convergentiecriterim is tweeledig. In de eerste plaats stopt de procedure, als er na een lijnminimalisatie geen

'±+1

f(x.)

uitkomt.

Na meer dan n lijnminimalisaties wordt bovendien op convergentie getest aan de hand van de volgende twee eisen:

(i) De linstaplengte x. - moet kleiner dan ziin. (2) De lengte van de zoekvector p. moet ook kleiner dan 10

zijn.

Als "lengte" of norm van een vector x neemt men hierbij

Il.H

=Max x.j.

(i)

De tweede eis betekent practisch zoveel als

Max g. < 10

-(i)

Dit dubbele criterium is nogal star.

Een gebruiker zal in het algemeen andere eisen stellen nl.

Hxi+1

-

-.x.II < E-i _x n f (x)jj < E

x

_g

Door enige modificaties In de oorspronkelijke functieprocedure aan te brengen, kan men zorgen, dat de, met optfpn geminimaliseerde functie

een minimuinpunt bereikt, dat aan de eisen, vervat

in

(149.1), voldoet.

Als we in de oorsprorikelijke functie f(x) voor x = O6 y inirullen, dan bet ekent

1j+1 -

< 1O, dat

-x.H

-i

<E

X (145.2) Het startpunt moeten we dan natuurlijk vervangen door

-= _6 Nu Is f(x) = f(106E ) f*() dus f (x) f*()1O_6/ X

af

6

a en =10 /(a

)-ay. ay. 3x. X ay.

i i i

< _{is, dan moet} af

x.

i

Dus

1O2/(ac

) = a =

Samenvatting:

Als men, gebruik makend van optfpn, wenst dat de procedure afbreekt op basis van het afbreekcriterium

x. -

xjj <c

n lf (x)J< E

-i+1 -J- x

x

g

dan moet men de oorspronkelike Í'unctieprocedure als voigt wizigen:

(i) Vervang x. door 10 c x.. (2) Vervang f() door

Vervang het startpunt door 10-6

Voorbeeld: Oorsoronkeiijke procedure: 'procedure'funct(x,fx);'array'x;'reai'fx; 'begin' 'real'xi ,x2;xl :=x(/1/) ;x2:=x(/2/); fx:=lOO*(xl-.x2*x2)**2+( 1-xi )**2 'funct; Gemodificeerd: 'procedure'funct(x,fx);'array'x;'real'fx; 'begin' 'reai'xfac,a,xl ,x2;xfac:epsx*'6;

a:=i/(epsx*ersg* '12); xl :=xfac*x(/1/) ;x2:=xfac*x(/2/); fx:a*(100*(xl_x2*x2)**2+(1_xl)**2)

''funct;

E

Zijfl.

g Als ay.

3. Minimalisering

van

functies met beperkingen.

3.1. Algemene probleeinstelling.

a is een willekeurige vector s Enì. _Op Rn als domein zijn gegeven de

reelwaardige functies f, h1 t/m h0 en t/m

De functies h.} en {g.}

zijn

zodanig, dat de verzameling punten B = {x:h.(x)=O,

i=1(1)p, g.(x)

O, j1(1)q}

niet leeg is. B is het z.g. realiseringsgebied.

Gevraagd wordt, het 'minimumpunt

van

f sub B17 te bepalen, d.w.z. het

punt x, waarvoor moet gelden dat

(i)

X E B

(ii) f() f(x) voor alle X B.

De eisen "gj(x) O" en "h.(x) = O" heten (on)gelijkheidsbeperkingen.

(Eng: (in)eq»uality constraints)

3 .2. Gelijkheidsbeperkingen.

Aan de hand van een eenvoudig voorbeeld zal worden verklaard aan welke

6o.rwaarden een punt x moet voldoen, dat minimumpunt is van f op de (n-p )-dimensionale deeiruimte

p

B= x

: h.2(x) = O}

i=1

Beschouw f(x) = X12

+ x1x2 +

en stel als enige eis

h(x)

= x1 +

-

2 = O.

1-We moeten in B - hier de lin h(x) f stationair is. Dit betekent, dat

een niveaulijn - aan B rnoet raken.

de gradient van f in x ioodrecht op B : h(x) = O.

In elk punt x staat dus h (x) loodrecht o B. De vectoren f (x) en moeten dus dezelfdeThrager hebben d.w.z. ze zijn iineir afhan-kelik. Er moet een getal zijn, zodanig dat

- T h(x) _{= Q.} In ons voorbeeld is f(.) = (2x1 + 2x2)t,

terwijl

h(.)

= (i i)'. Nu is dus 2x1 + = + 2 = , terwiji h(x) = O betekent, dat = = 2. Hieruit voigt = (1

1)', u

= 3. Lagrange-functie:

Men kan oDrnerken, dat deze opiossing ook verkregen kan worden, als men de z.g. Lagrange-functie

F(x,p) = f(.) - p h(x)

ininimaliseert t.a.v.

X en inaximaliseert t.a.v. p.

Uit

F( .) =

O voigt nl.

2x1

+ X2 -

p 0 en

x1

-

=

12

zodat F(x,p)

= -

i-t + 2p.

Deze functie van p is maximnal voor p = 3.

Merk op dat de eis "F(.) = O betekent, dat h(.) = OEmoet zijn. Bij ineerdere geiijkheidsrestricties

h(x)

O, i

= 1(1)p,

=

O -

een punt opsporen, waarin een niveau-opperviak van

f -

hier

Dit "raken aan B" betekent dat

B moet staan.

Nu geldt voor

moet in

het ininimumpunt x van

f sub

B voidaan zijn

aan

)49

-d.w.z. de gradint van f moet lineair afhankelijk zijn van de gradiënten van h1 t/m h

Voert men in:

h = (h1 . h)' en i = ...

dan is de eis (5.3.2) oak te schrijven als

f () - h ()

X-

_-X

-

=

-(n,1)

_{(n,p) (p,i)}

Hier kan warden opgemerkt dat het minimumpunt x en de bi.jbehorende

Lagrange-vector ii oak kunnen warden gevonden door de Lagrange-functie

F(x,) = f(x) - _(19.2)

ten aanzien van x te minimaliseren en t.a.v. te maximaliseren.

Een voldoande voorwaarde voor een minimumpunt van F is, dat de Hesse-matrix an F positief definiet is. Deze Hesse-matrix is hier een (n,n)-matrix met in de e rij en e kolom het element

2f -X_ x. x. i j k=1 i j

(91)

(19.3)