Stochastické algoritmy pro globální optimalizaci

Stochastick´ e algoritmy pro glob´ aln´ı optimalizaci Josef Tvrd´ık

ˇ ÍSLO OPERACN ˇ ÍHO PROGRAMU: CZ.1.07 C ´ ˇ ÍHO PROGRAMU: NAZEV OPERACN ˇ AV ´ AN ´ Í PRO KONKURENCESCHOPNOST OP VZDEL PRIORITNÍ OSA: 2 ˇ ÍSLO OBLASTI PODPORY: 2.3 C POSÍLENÍ KONKURENCESCHOPNOSTI ´ ´ ˇ ÍCH TECHNOLOGIÍ VYZKUMU A VYVOJE INFORMACN ´ V MORAVSKOSLEZSKEM KRAJI ˇ Í C ˇ ÍSLO PROJEKTU: CZ.1.07/2.3.00/09.0197 REGISTRACN

OSTRAVA 2010

Tento projekt je spolufinancován Evropsk´ ym sociáln´ım fondem a státn´ım rozpoˇctem ˇ e republiky Cesk´

Název: Autor:

Stochastické algoritmy pro globáln´ı optimalizaci Josef Tvrd´ık

Vydán´ı: Poˇcet stran:

prvn´ı, 2010 80

Studijn´ı materiály pro distanˇcn´ı kurz: Stochastické algoritmy pro globáln´ı optimalizaci

Jazyková korektura nebyla provedena, za jazykovou stránku odpov´ıdá autor.

c Josef Tvrd´ık ⃝ c Ostravská univerzita v Ostravˇe ⃝

OBSAH

1

Obsah ´ 1 Uvod

2

2 Probl´ em glob´ aln´ı optimalizace

3

2.1

Formulace problému globáln´ı optimalizace . . . . . . . . . . . . . . .

3

2.2

Stochastické algoritmy pro globáln´ı optimalizaci . . . . . . . . . . . .

6

3 Experiment´ aln´ı porovn´ av´ an´ı algoritm˚ u

7

3.1

Experimentáln´ı ovˇeˇrován´ı algoritm˚ u. . . . . . . . . . . . . . . . . . .

7

3.2

Testovac´ı funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4 Slep´ e n´ ahodn´ e prohled´ av´ an´ı

20

ˇ ızen´ 5 R´ e n´ ahodn´ e prohled´ av´ an´ı

23

6 Evoluˇ cn´ı algoritmy

30

6.1

Genetické algoritmy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

6.2

Evoluˇcn´ı strategie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

6.3

Diferenciáln´ı evoluce . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

6.4

Experimentáln´ı porovnán´ı jednoduch´ ych evoluˇcn´ıch algoritm˚ u . . . . 52

7 Algoritmy modeluj´ıc´ı chov´ an´ı skupin

55

7.1

Algoritmus PSO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

7.2

Algoritmus SOMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

8 Adaptace pomoc´ı soutˇ eˇ ze strategi´ı

70

9 Literatura

78

´ 1 UVOD

2

1

´ Uvod

Tento text je urˇcen student˚ um volitelného pˇredmˇetu Stochastické algoritmy pro globáln´ı optimalizaci v prezenˇcn´ım, kombinovaném i v distanˇcn´ım studiu na Ostravské univerzitˇe. C´ılem textu je seznámit studenty se základy problematiky globáln´ı optimalizace a algoritmy pro heuristické hledán´ı globáln´ıho extrému funkc´ı. Text je zamˇeˇren pˇredevˇs´ım na stochastické algoritmy inspirované ˇzivou pˇr´ırodou, tj. na evoluˇcn´ı algoritmy a na algoritmy modeluj´ıc´ı chován´ı spoleˇcenstv´ı ˇzivoˇciˇsn´ ych jedinc˚ u v kolektivu. Kaˇzdá kapitola zaˇc´ıná pokyny pro jej´ı studium. Tato ˇca´st je vˇzdy oznaˇcena jako Pr˚ uvodce studiem s ikonou na okraji stránky. Pojmy a d˚ uleˇzité souvislosti k zapamatován´ı jsou vyznaˇceny na okraji stránky textu ikonou. V závˇeru kaˇzdé kapitoly je rekapitulace nejd˚ uleˇzitˇejˇs´ıch pojm˚ u. Tato rekapitulace je oznaˇcena textem Shrnut´ı a ikonou na okraji. Odd´ıl Kontroln´ı ot´ azky oznaˇcen´ y ikonou by vám mˇel pomoci zjistit, zda jste prostudovanou kapitolu pochopili a budou i inspirovat k dalˇs´ım otázkám, na které budete hledat odpovˇed’. U nˇekter´ ych kapitol je pˇripomenuto, ˇze k této ˇcásti textu je zadána Korespon´ eˇsné vyˇreˇsen´ı korespondenˇcn´ıch u denˇ cn´ı u ´ loha. Uspˇ ´loh je souˇca´st´ı podm´ınek pro z´ıskán´ı zápoˇctu pro studenty kombinovaného a distanˇcn´ıho studia. Korespondenˇcn´ı u ´lohy budou zadávány v rámci kurzu daného semestru. V textu budeme uˇz´ıvat následuj´ıc´ı druhy p´ısma pro r˚ uzné symboly: • kurz´ıva malá p´ısmena, napˇr. a, b, x1 , yj , δ, ε pro re´ aln´ a ˇc´ısla, • kurz´ıva velká p´ısmena, napˇr. X, S, D pro mnoˇziny reáln´ ych ˇc´ısel nebo vektor˚ u. Pro celou mnoˇzinu reáln´ ych ˇc´ısel uˇz´ıváme speciáln´ı symbol R, • kurz´ıva tuˇcná malá p´ısmena, napˇr. a, b, x, α pro vektory. Vektor délky n je uspoˇrádaná n-tice reáln´ ych ˇc´ısel, napˇr. x = (x1 , x2 , . . . , xn ), • kurz´ıva tuˇcná velká p´ısmena, napˇr. X, S pro matice. Zdrojové texty algoritm˚ u v Matlabu jsou tiˇstˇeny strojopisn´ ym typem, kter´ y vypadá napˇr. takto: y = x + (b-a) .* rand(1,n). Pokud nejste obeznámeni se základy práce s Matlabem, doporuˇcujeme nejdˇr´ıve absolvovat kurz Základy modelován´ı v MATLABU nebo alespoˇ n prostudovat uˇcebn´ı texty k tomuto kurzu [24].

3

2

Probl´ em glob´ aln´ı optimalizace

Pr˚ uvodce studiem: V této kapitole je zformulován problém globáln´ı optimalizace a struˇcnˇe jsou zm´ınˇeny základn´ı myˇslenky stochastick´ ych algoritm˚ u, které se pouˇz´ıvaj´ı pro nalezen´ı globáln´ıho minima u ´ˇcelov´ ych funkc´ı. Poˇc´ıtejte asi se dvˇema hodinami studia.

2.1

Formulace probl´ emu glob´ aln´ı optimalizace

Budeme se zab´ yvat ˇreˇsen´ım problému globáln´ı optimalizace, tj. nalezen´ım souˇradnic takového bodu v definiˇcn´ım oboru funkce, ve kterém má funkce extrémn´ı (nejmenˇs´ı nebo nejvˇetˇs´ı) hodnotu. Problém nalezen´ı globáln´ıho minima funkce s jedn´ım argumentem (argument je jedno reálné ˇc´ıslo) m˚ uˇzeme ilustrovat jednoduch´ ym obrázkem:

1

f(x)

0.5

0

−0.5 globální minimum x* −1 0

5

10 D=[0; 18], d = 1

15

Vid´ıme, ˇze funkce na obrázku má v oboru [0, 18] v´ıce minim, ale jen jedno je globáln´ı, tj. takové, ˇze funkˇcn´ı hodnota v tomto bodˇe je nejmenˇs´ı ze vˇsech hodnot v oboru [0, 18]. Ze základn´ıho kurzu matematické anal´ yzy známe postup, jak nalézt extrémy funkc´ı, u kter´ ych existuje prvn´ı a druhá derivace. Zdálo by se, ˇze u ´loha nalezen´ı globáln´ıho minima je velmi jednoduchá. Bohuˇzel tomu tak nen´ı. Nalézt obecné ˇreˇsen´ı takto jednoduˇse pochopitelného problému je obt´ıˇzné, zvláˇstˇe kdyˇz u ´ˇcelová funkce má v´ıce lokáln´ıch minim nebo argument funkce nen´ı jedno reálné ˇc´ıslo, ale vektor reáln´ ych ˇc´ısel, viz dalˇs´ı obrázek pro funkci se dvˇema argumenty. Nav´ıc, ne kaˇzdá funkce

´ GLOBALN ´ Í OPTIMALIZACE 2 PROBLEM

4

je diferencovatelná, ale pˇresto potˇrebujeme jej´ı globáln´ı extrém nalézt nebo se mu alespoˇ n pˇribl´ıˇzit s pˇrijatelnou pˇresnost´ı.

1000

f(x)

500

0

−500

−1000 −500 −500 0 0 500

500

x2

Global minimum point x1

´ Ulohu nalezen´ı globáln´ıho minima m˚ uˇzeme formulovat takto: Mˇejme u ´ˇcelovou funkci f : D → R,

D ⊆ Rd .

(2.1)

Máme naj´ıt bod x∗ ∈ D, pro kter´ y plat´ı, ˇze f (x∗ ) ≤ f (x), pro ∀x, x ∈ D. Nalezen´ı bodu x∗ ∈ D je ˇreˇsen´ım problému globáln´ı optimalizace. Bodu x∗ ˇr´ıkáme bod globáln´ıho minima (global minimum point), definiˇcn´ımu oboru D se ˇr´ıká doména nebo prohledávan´ y prostor (domain, search space), pˇrirozené ˇc´ıslo d je dimenze u ´lohy. Formulace problému globáln´ı optimalizace jako nalezen´ı globáln´ıho minima nen´ı na u ´kor obecnosti, nebot’ chceme-li nalézt globáln´ı maximum, pak jej nalezneme jako globáln´ı minimum funkce g(x) = −f (x). Anal´ yza problému globáln´ı optimalizace ukazuje, ˇze neexistuje deterministick´ y algoritmus ˇreˇs´ıc´ı obecnou u ´lohu globáln´ı optimalizace (tj. nalezen´ı dostateˇcnˇe pˇresné ∗ aproximace x ) v polynomiáln´ım ˇcase, tzn. problém globáln´ı optimalizace je NPobt´ıˇzn´ y. Pˇritom globáln´ı optimalizace je u ´loha, kterou je nutno ˇreˇsit v mnoha praktick´ ych problémech, mnohdy s velmi v´ yznamn´ ym ekonomick´ ym efektem, takˇze je nutné hledat algoritmy, které jsou pro ˇreˇsen´ı konkrétn´ıch problém˚ u pouˇzitelné. Algoritmy pro ˇreˇsen´ı problému globáln´ı optimalizace se podrobnˇe zab´ yvá celá ˇrada monografi´ı, ˇ ych napˇr. Torn a Zilinskas [22], M´ıka [15], Spall [19], kde je moˇzné naj´ıt mnoho uˇziteˇcn´ poznatk˚ u pˇresahuj´ıc´ıch rámec tohoto pˇredmˇetu.

2.1 Formulace problému globáln´ı optimalizace

5

´ Uloha (2.1) se oznaˇcuje jako hledán´ı volného extrému funkce (unconstraint optimization). Je moˇzné pˇrijatelnost ˇreˇsen´ı jeˇstˇe omezit (vázat) nˇejakou podm´ınkou, napˇr. nˇejak´ ymi rovnicemi nebo nerovnostmi. Pak jde o problém hledán´ı vázaného extrému (constrained optimization). V tomto kurzu se soustˇred´ıme pˇredevˇs´ım na problémy, kdy hledáme globáln´ı minimum v souvislé oblasti D = ⟨a1 , b1 ⟩ × ⟨a2 , b2 ⟩ × . . . × ⟨ad , bd ⟩ =

∏d

i=1 ⟨ai , bi ⟩,

(2.2)

ai < bi , i = 1, 2, . . . , d, au ´ˇcelovou funkci f (x) um´ıme vyhodnotit s poˇzadovanou pˇresnost´ı v kaˇzdém bodˇe x ∈ D. Podm´ınce (2.2) se ˇr´ıká boundary constraints nebo box constraints, protoˇze oblast D je vymezena jako d-rozmˇern´ y kvádr. Pro u ´lohy ˇreˇsené numericky na poˇc´ıtaˇci nepˇredstavuje podm´ınka (2.2) ˇza´dné podstatné omezen´ı, nebot’ hodnoty ai , bi jsou tak jako tak omezeny datov´ ymi typy uˇzit´ ymi pro x a f (x), tj. vˇetˇsinou reprezentac´ı ˇc´ısel v pohyblivé ˇrádové ˇca´rce. Proto se dosti ˇcasto takové u ´lohy oznaˇcuj´ı jako unconstrained continuous problems. ´ Ulohy hledán´ı vázaného extrému (constrained optimization) jsou obvykle formulovány takto: Najdi minimum funkce f (x), x = (x1 , x2 , . . . , xd ) a x ∈ D

(2.3)

za podm´ınky: gi (x) ≤ 0, i = 1, . . . , p hj (x) = 0, j = p + 1, . . . , m. ˇ sen´ı je povaˇzováno za pˇrijatelné (feasible), kdyˇz gi (x) ≤ 0 pro i = 1, . . . , p a Reˇ |hj (x)| − ε ≤ 0 pro j = p + 1, . . . , m. Pro libovoln´ y bod x ∈ D a zadané kladné ˇc´ıslo ε m˚ uˇzeme definovat pr˚ umˇerné poruˇsen´ı podm´ınek (mean violation) v¯ jako ∑p v¯ = kde

∑m G (x) + i j=p+1 Hj (x) i=1 , m

  gi (x) Gi (x) =  0   |hj (x)| Hj (x) =  0

if

gi (x) > 0

if

gi (x) ≤ 0

if

|hj (x)| − ε > 0

if

|hj (x)| − ε ≤ 0.

Moˇzná jste se s nejjednoduˇsˇs´ı variantou takov´ ych u ´loh setkali. Omezuj´ıc´ı podm´ınkou byla lineárn´ı nerovnost a u ´lohu jste vyˇreˇsili tzv. lineárn´ım programován´ım.

6

´ GLOBALN ´ Í OPTIMALIZACE 2 PROBLEM

Dalˇs´ım typem optimalizaˇcn´ıch u ´loh jsou diskrétn´ı problémy, kdy prohledáváná oblast D nen´ı spojitá, ale diskrétn´ı, napˇr. hodnoty jednotliv´ ych prvk˚ u vektoru x mohou b´ yt pouze celoˇc´ıselné. Do této skupiny u ´loh patˇr´ı napˇr. problémy hledán´ı optimáln´ı cesty v grafu.

2.2

Stochastick´ e algoritmy pro glob´ aln´ı optimalizaci

Nemoˇznost nalézt deterministick´ y algoritmus obecnˇe ˇreˇs´ıc´ı u ´lohu globáln´ı optimalizace v polynomiáln´ım ˇcase vedla k vyuˇzit´ı algoritm˚ u stochastick´ ych, které sice nemohou garantovat nalezen´ı ˇreˇsen´ı v koneˇcném poˇctu krok˚ u, ale ˇcasto pomohou nalézt v pˇrijatelném ˇcase ˇreˇsen´ı prakticky pouˇzitelné. Stochastické algoritmy pro globáln´ı optimalizaci heuristicky prohledávaj´ı prostor D. Heuristikou rozum´ıme postup, ve kterém se vyuˇz´ıvá náhoda, intuice, analogie a zkuˇsenost. Rozd´ıl mezi heuristikou a deterministick´ ym algoritmem je v tom, ˇze na rozd´ıl od deterministického algoritmu heuristika nezajiˇst’uje nalezen´ı ˇreˇsen´ı. Heuristiky jsou v praktickém ˇzivotˇe zcela samozˇrejmˇe uˇz´ıvané postupy, jako pˇr´ıklady m˚ uˇzeme uvést hledán´ı hub, lov ryb na udici, v´ ybˇer partnera, pokus o v´ yhru ve sportovn´ım utkán´ı nebo o sloˇzen´ı zkouˇsky ve ˇskole. Vˇetˇsina stochastick´ ych algoritm˚ u pro hledán´ı globáln´ıho minima v sobˇe obsahuje zjevnˇe ˇci skrytˇe proces uˇcen´ı. Inspirace k uˇzit´ı heuristik jsou ˇcasto odvozeny ze znalost´ı pˇr´ırodn´ıch nebo sociáln´ıch proces˚ u. Napˇr. simulované ˇz´ıhán´ı je modelem pomalého ochlazován´ı tuhého tˇelesa, tabu-search modeluje hledán´ı pˇredmˇetu ˇclovˇekem tak, ˇze v krátkodobé pamˇeti si zapamatovává zakázané kroky vedouc´ı k jiˇz dˇr´ıve projit´ ym m´ıst˚ um. Popis tˇechto algoritm˚ u najdete napˇr. v knize Kvasniˇcky, Posp´ıchala a Tiˇ na [12]. Podobné postupy uˇcen´ı najdeme snad ve vˇsech znám´ ych stochastick´ ych algoritmech pro hledán´ı globáln´ıho minima (s v´ yjimkou slepého náhodného prohledáván´ı). Znaˇcná ˇcást stochastick´ ych algoritm˚ u pracuje souˇcasnˇe s v´ıce kandidáty ˇreˇsen´ı, tj. s v´ıce body v prohledávaném prostoru. Tyto body vytváˇrej´ı skupinu (populaci), která se v pr˚ ubˇehu hledán´ı nˇejak v prohledávaném prostoru pohybuje a pˇritom nacház´ı lepˇs´ı kandidáty ˇreˇsen´ı. Stochastické algoritmy pro globáln´ı optimalizaci jsou pˇr´ıkladem vyuˇzit´ı soft computingu pro ˇreˇsen´ı u ´loh, které hard computingem ˇreˇsit neum´ıme.

7

3

Experiment´ aln´ı porovn´ av´ an´ı algoritm˚ u

Pr˚ uvodce studiem: V této kapitole je uvedeno nˇekolik funkc´ı, které se uˇz´ıvaj´ı k testován´ı stochastick´ ych algoritm˚ u pro globáln´ı optimalizaci. Z nich si vybereme testovac´ı funkce pro experimenty. Druhá ˇca´st kapitoly se zab´ yvá testován´ım stochastick´ ych algoritm˚ u. K orientaci v této kapitole budete potˇrebovat asi dvˇe aˇz tˇri hodiny.

3.1

Experiment´ aln´ı ovˇ eˇ rov´ an´ı algoritm˚ u

Pro ovˇeˇrován´ı a vzájemné porovnáván´ı stochastick´ ych algoritm˚ u jsou nutné numerické experimenty. Pˇri porovnáván´ı dvou ˇci v´ıce algoritm˚ u je samozˇrejmˇe nutné, aby testy byly provedeny za stejn´ ych podm´ınek, tj. zejména pˇri stejnˇe vymezeném prohledávaném prostoru D a pˇri stejné podm´ınce pro ukonˇcen´ı prohledáván´ı. Základn´ımi veliˇcinami pro porovnán´ı algoritm˚ u jsou ˇcasová nároˇcnost prohledáván´ı a ˇ spolehlivost nalezen´ı globáln´ıho minima. Casov´ a nároˇcnost algoritmu se ohodnocuje poˇctem vyhodnocen´ı u ´ˇcelové funkce v pr˚ ubˇehu hledán´ı, takˇze v´ ysledky jsou srovnatelné bez ohledu na rychlost poˇc´ıtaˇc˚ u uˇzit´ ych k testován´ı. Poˇcet vyhodnocen´ı u ´ˇcelové funkce potˇrebn´ y k dosaˇzen´ı podm´ınky ukonˇcen´ı m˚ uˇzeme oznaˇcit napˇr. nfe (number of function evaluations). U testovac´ıch funkc´ı, kdy ˇreˇsen´ı problému je známé, se nˇekdy také sleduje poˇcet vyhodnocen´ı u ´ˇcelové funkce nfe_near potˇrebn´ y k tomu, aby nejlepˇs´ı bod populace (skupiny) mˇel hodnotu niˇzˇs´ı neˇz je zadaná hodnota fnear , o které se v´ı, ˇze jej´ı dosaˇzen´ı znamená, ˇze se algoritmus nejlepˇs´ım nalezen´ ym bodem dostateˇcnˇe pˇribl´ıˇzil globáln´ımu minimu. Nˇekdy se tato hodnota naz´ yvá také VTR, coˇz je zkratka anglického term´ınu value to reach. Podm´ınka ukonˇcen´ı b´ yvá formulována jako fmax − fmin < ε, kde fmax je nejvˇetˇs´ı a fmin nejmenˇs´ı funkˇcn´ı hodnota mezi vˇsemi body skupiny (populace). Prohledáván´ı konˇc´ı, kdyˇz rozd´ıl funkˇcn´ıch hodnot je tak mal´ y, ˇze uˇz nelze oˇcekávat nalezen´ı vhodnˇejˇs´ıch kandidát˚ u ˇreˇsen´ı. Je vˇsak praktické formulovat podm´ınku ukonˇcen´ı tak, abychom pˇredeˇsli abnormáln´ımu ukonˇcen´ı programu napˇr. v situaci, kdy prohledáván´ı dojde do nekoneˇcného cyklu. Proto je vhodné pˇridat dalˇs´ı vstupn´ı parametr algoritmu – maxevals, coˇz

8

´ Í POROVNAV ´ AN ´ Í ALGORITM˚ 3 EXPERIMENTALN U

je maximáln´ı dovolen´ y poˇcet vyhodnocen´ı funkce na jednu dimenzi prohledávaného prostoru a podm´ınku ukonˇcen´ı pak formulovat jako fmax − fmin < ε ∨ nfe ≥ maxevals ∗ d, to znamená, ˇze hledán´ı pokraˇcuje tak dlouho, dokud funkˇcn´ı hodnoty v populaci se liˇs´ı v´ıce neˇz poˇzadujeme nebo dokud nen´ı dosaˇzen nejvyˇsˇs´ı dovolen´ y poˇcet vyhodnocen´ı u ´ˇcelové funkce, tzn. je splnˇena podm´ınka fmax − fmin ≥ ε ∧ nfe < maxevals ∗ d. M´ısto fmax lze v podm´ınce ukonˇcen´ı uˇz´ıt i jin´ y kvantil funkˇcn´ı hodnoty v populaci, napˇr. medián. T´ım uˇcin´ıme podm´ınku ukonˇcen´ı snadnˇeji dosaˇzitelnou, staˇc´ı, aby jen napˇr. polovina populace mˇela funkˇcn´ı hodnoty dostateˇcnˇe bl´ızké. Z uveden´ ych u ´vah je zˇrejmé, ˇze pak mus´ıme rozliˇsovat následuj´ıc´ı ˇctyˇri typy ukonˇcen´ı prohledáván´ı v tetovac´ıch u ´lohách: • Typ 1 – korektn´ı ukonˇcen´ı, tj. fmax − fmin < ε a fmin ≤ fnear , algoritmus tedy splnil podm´ınku ukonˇcen´ı pˇred dosaˇzen´ım maximáln´ıho dovoleného poˇctu iterac´ı a souˇcasnˇe se dostateˇcnˇe pˇribl´ıˇzil globáln´ımu minimu. • Typ 2 – pomalá konvergence, fmin ≤ fnear , ale prohledáván´ı bylo ukonˇceno dosaˇzen´ım maximáln´ıho dovoleného poˇctu iterac´ı. • Typ 3 – pˇredˇcasné konvergence (premature convergence), fmax − fmin < ε, ale fmin > fnear , tzn. nebylo nalezeno globáln´ı minimum, algoritmus bud’ skonˇcil prohledáván´ı v lokáln´ım minimu nebo se body populace pˇribl´ıˇzily k sobˇe natolik, ˇze jejich funkˇcn´ı hodnoty jsou velmi bl´ızké. • Typ 4 – u ´plné selhán´ı, prohledáván´ı bylo ukonˇceno dosaˇzen´ım maximáln´ıho dovoleného poˇctu iterac´ı, aniˇz byl nalezen bod bl´ızk´ y globáln´ımu minimu. Po prohledávac´ıch algoritmech globáln´ı optimalizace chceme, aby umˇely nacházet ˇreˇsen´ı dostateˇcnˇe bl´ızké globáln´ımu minimu x∗ co nejrychleji, co nejspolehlivˇeji a aby prohledáván´ı umˇely ukonˇcit ve vhodném stádiu prohledáván´ı. Zcela u ´spˇeˇsné prohledáván´ı je tedy jen takové, jehoˇz ukonˇcen´ı je typu 1, nˇekdy vˇsak m˚ uˇzeme povaˇzovat za u ´spˇeˇsné i ukonˇcen´ı typu 2. Spolehlivost (reliability, R) nalezen´ı globáln´ıho minima m˚ uˇzeme charakterizovat jako relativn´ı ˇcetnost ukonˇcen´ı typu 1 (pˇr´ıpadnˇe typu 1 nebo typu 2) R=

n1 , n

(3.1)

kde n1 je poˇcet ukonˇcen´ı typu 1 (pˇr´ıpadnˇe typu 1 nebo 2) v n nezávisl´ ych opakován´ıch. Veliˇcina n1 je náhodná veliˇcina n1 ∼ Bi(n, p) a R je nestrann´ ym odhadem

3.1 Experimentáln´ı ovˇeˇrován´ı algoritm˚ u

9

parametru p, tj. pravdˇepodobnosti, ˇze algoritmus v daném problému najde globáln´ı minimum (pˇresnˇeji, ˇze konec prohledáván´ı je typu 1). Rozptyl var(R) = p(1 − p)/n. Pro( vˇetˇs´ı hodnoty n m˚ uˇzeme rozdˇelen´ı R povaˇzovat za pˇribliˇznˇe normáln´ı, R ∼ ) p(1−p) N p, n . Pak 100 (1 − α)-procentn´ı oboustrann´ y interval spolehlivosti pro p je [R − △, R + △], √ α) R(1 − R) △ = tn−1 1 − 2 n V následuj´ıc´ı tabulce jsou uvedeny hodnoty △ pro α = 0.05 a pro r˚ uzná n a R. Pokud byste vyjadˇrovali spolehlivost R v procentech, pak je potˇreba hodnoty △ v tabulce vynásobit stovkou. Z uvedené tabulky je zˇrejmé, ˇze pro porovnáván´ı spolehlivosti algoritm˚ u potˇrebujeme velk´ y poˇcet opakován´ı n. Napˇr. pro n = 100 je moˇzno za v´ yznamn´ y povaˇzovat aˇz rozd´ıl v reliabilitˇe vˇetˇs´ı neˇz 10 aˇz 20 procent (závis´ı na hodnotˇe R). kde

(

Hodnoty △ pro α = 0.05 R n

0.50

0.60

0.70

0.80

0.90

0.95

10 0.358 0.350 0.328 0.286 0.215 0.156 25 0.206 0.202 0.189 0.165 0.124 0.090 50 0.142 0.139 0.130 0.114 0.085 0.062 100 0.099 0.097 0.091 0.079 0.060 0.043 200 0.070 0.068 0.064 0.056 0.042 0.030 500 0.044 0.043 0.040 0.035 0.026 0.019 V´ ysledky testován´ı by vˇzdy mˇely obsahovat specifikaci testovan´ ych algoritm˚ u vˇcetnˇe hodnot jejich vstupn´ıch parametr˚ u, specifikaci testovac´ıch funkc´ı, poˇcet opakován´ı experiment˚ u, tabulku s pr˚ umˇern´ ymi hodnotami charakterizuj´ıc´ımi ˇcasovou nároˇcnost, pˇr´ıpadnˇe i vhodné charakteristiky jej´ı variability a spolehlivost R, která m˚ uˇze b´ yt uvedena v procentech. Pro názorné porovnán´ı ˇcasové nároˇcnosti jsou vhodné krabicové grafy, ve kter´ ych snadno vizuálnˇe porovnáme jak charakteristiky polohy, tak variability. Pokud rozd´ıly v charakteristikách algoritm˚ u nejsou zjevné z popisné statistiky, pouˇzijte vhodné statistické testy. Pro porovnán´ı ˇcasové nároˇcnosti to m˚ uˇze b´ yt dvouv´ ybˇerov´ y t -test, anal´ yza rozptylu nebo jejich neparametrické analogie. Pokud chcete porovnávat spolehlivosti dvou algoritm˚ u, m˚ uˇzete pro ˇctyˇrpoln´ı tabulku absolutn´ıch ˇcetnost´ı uˇz´ıt Fisher˚ uv exaktn´ı test.

10


U nˇekter´ ych algoritm˚ u je zaj´ımavé sledovat i jiné veliˇciny, které charakterizuj´ı pr˚ ubˇeh prohledávac´ıho procesu. Volba dalˇs´ıch sledovan´ ych veliˇcin je závislá na testovaném algoritmu. Zaj´ımavou popisnou informac´ı o pr˚ ubˇehu vyhledáván´ı je také graf závislosti funkˇcn´ıch hodnot na poˇctu vyhodnocen´ı u ´ˇcelové funkce. Pr˚ ubˇeh této závislosti nám umoˇzn ˇuje posoudit rychlost konvergence algoritmu v r˚ uzn´ ych fáz´ıch prohledáván´ı. V takov´ ych grafech je ovˇsem ˇzádouc´ı m´ıt charakteristiky z v´ıce opakován´ı a pokud moˇzno nejen pr˚ umˇerné hodnoty, ale i minima a maxima. K rychlému porovnáván´ı u ´ˇcinnosti algoritm˚ u byly navrˇzeny r˚ uzné dalˇs´ı veliˇciny integruj´ıc´ı ˇcasovou nároˇcnost i spolehlivost do jediného kritéria. Jednou z nich je tzv. Q-m´ıra, definovaná jako Qm = nfe/R, kde nfe je pr˚ umˇern´ y poˇcet vyhodnocen´ı funkce v n opakovan´ ych bˇez´ıch algoritmu pˇri ˇreˇsen´ı u ´lohy a R je spolehlivost definovaná vztahem (3.1). Pak nej´ uˇcinnˇejˇs´ı algoritmus mezi porovnávan´ ymi je ten s minimáln´ı hodnotou Qm . Nev´ yhodou tohoto kritéria je, ˇze jeho hodnota nen´ı definována pro R = 0. Proto je vhodnˇejˇs´ı Qinv = R/nfe, pro kterou plat´ı, ˇze ˇc´ım vyˇsˇs´ı hodnota Qinv , t´ım lepˇs´ı u ´ˇcinnost algoritmu. Kritéria integruj´ıc´ı nfe a R jsou vˇsak vhodná jen pro rychlé pˇrehledné porovnáván´ı celkové u ´ˇcinnosti algoritm˚ u, protoˇze z nich nelze dˇelat u ´sudky o variabilitˇe ˇcasové nároˇcnosti algoritm˚ u, ani jednoznaˇcnˇe posoudit jejich spolehlivost.

3.2

Testovac´ı funkce

K testován´ı stochastick´ ych algoritm˚ u pro globáln´ı optimalizaci se uˇz´ıvá ˇrada funkc´ı, u kter´ ych je známé správné ˇreˇsen´ı, tj. globáln´ı minimum, které lze naj´ıt analyticky. Testovac´ı funkce (benchmark functions) jsou ˇcasto navrˇzeny tak, aby je bylo moˇzné pouˇz´ıt pro problémy s r˚ uznou dimenz´ı prohledávané domény. Testovac´ı funkce mohou m´ıt r˚ uznou obt´ıˇznost, nejlehˇc´ı jsou konvexn´ı funkce, ty maj´ı jen jedno lokáln´ı minimum, které je tedy i minimem globáln´ım. Takov´ ym funkc´ım se ˇr´ıká unimodáln´ı. Funkce s v´ıce lokáln´ımi minimy se naz´ yvaj´ı multimodáln´ı a nalezen´ı globáln´ıho minima je obvykle obt´ıˇznˇejˇs´ı neˇz u unimodáln´ıch funkc´ı. Za lehˇc´ı problémy jsou povaˇzována hledán´ı globáln´ıho minima separabiln´ıch funkc´ı, tj. funkc´ı, jejichˇz argumenty jsou nezávislé.

3.2 Testovac´ı funkce

11

Obecnˇe plat´ı, ˇze funkce splˇ nuj´ıc´ı podm´ınku ∂f (x) = g(xi ) h(x) ∂xi jsou separabiln´ı. Pro takové funkce plat´ı ∀i, j

i ̸= j :

f (. . . , x∗i , . . .) = opt ∧ f (. . . , x∗j , . . .) = opt ⇒ f (. . . , x∗i , . . . , x∗j , . . .) = opt

ˇ Casto lze takové funkce vyhodnotit jako f ((x) =

d ∑

fi (xi )

(3.2)

i=1

Hledán´ı minima takov´ ych separabiln´ıch funkc´ı je snazˇs´ı, protoˇze poloha x∗i je nezávislá na poloze x∗j pro i ̸= j a tedy mohou b´ yt hledány nezávisle na sobˇe. Zde uvád´ıme nˇekolik funkc´ı (separabiln´ıch i neseparabiln´ıch), které m˚ uˇzeme uˇz´ıt pˇri experimentáln´ı testován´ı algoritm˚ u. Dalˇs´ı testovac´ı funkce naleznete napˇr. v [2, 20, 29] nebo na webu. Vˇsechny dále uvedené funkce lze uˇz´ıt pro r˚ uzné dimenze prohledávané domény pro d ≥ 2, nˇekteré i pro d = 1.


12

Prvn´ı De Jongova funkce je paraboloid f (x) =

d ∑

x2i ,

i=1

pˇri experimentech prohledávan´ y prostor je obvykle vymezen podm´ınkou xi ∈ [−5.12, 5.12]. Je to konvexn´ı separabiln´ı funkce, globáln´ı minimum je x∗ = (0, . . . , 0), f (x∗ ) = 0. Je povaˇzována za velmi snadnou u ´lohu globáln´ı optimalizace. Pokud algoritmus selhává nebo pomalu konverguje v této u ´loze, tak nezaslouˇz´ı dalˇs´ı pozornosti. Graf této funkce pro d = 2 je na následuj´ıc´ım obrázku.

50

40

30

20

10

0 5 5 0 0 −5

−5


13

Druhá De Jongova funkce (Rosenbrockovo sedlo, známá i pod názvem banánové u ´dol´ı), neseparabiln´ı funkce, je definována takto: f (x) =

d−1 ∑ [ ] 100(x2i − xi+1 )2 + (1 − xi )2 , i=1

prohledávan´ y prostor je obvykle vymezen podm´ınkou xi ∈ [−2.048, 2.048], globáln´ı minimum je x∗ = (1, 1, . . . 1), f (x∗ ) = 0. Aˇc je to jednomodáln´ı funkce (v posledn´ı dobˇe se objevily domnˇenky, ˇze pro nˇekteré dimenze má tato funkce nev´ yrazná lokáln´ı minima), nalezen´ı minima iteraˇcn´ımi algoritmy nen´ı jednoduchá u ´loha, nebot’ minimum leˇz´ı v zahnutém u ´dol´ı s velmi mal´ ym spádem a nˇekteré algoritmy konˇc´ı prohledáván´ı pˇredˇcasnˇe. Graf této funkce pro d = 2 je na následuj´ıc´ım obrázku.

3000

2000 minimum f [1,1] = 0 1000

0 4 2 0 −2

−2

−1

0

1

2


14

Ackleyho funkce je v  ( d ) u d u1 ∑ ∑ 1 f (x) = −20 exp −0.02 t x2  − exp cos 2πxi + 20 + exp(1), d i=1 i d i=1 

xi ∈ [−30, 30], x∗ = (0, 0, . . . 0), f (x∗ ) = 0. Je to multimodáln´ı separabiln´ı funkce s mnoha lokáln´ımi minimy, stˇrednˇe obt´ıˇzná. Graf této funkce pro d = 2, xi ∈ [−3, 3], je na následuj´ıc´ım obrázku.

4 3

f(x)

2 1 0 −1 2 0 −2 x2

−3

−2

0

−1 x1

1

2

3


15

Griewankova funkce, multimodáln´ı, neseparabiln´ı, ( ) d d ∑ ∏ xi 2 xi f (x) = − cos √ + 1, 4000 i=1 i i=1 xi ∈ [−400, 400], x∗ = (0, 0, . . . 0), f (x∗ ) = 0. Nalezen´ı minima této multimodáln´ı funkce je povaˇzováno za obt´ıˇznou u ´lohu globáln´ı optimalizace. Graf této funkce pro d = 2, xi ∈ [0, 15], je na následuj´ıc´ım obrázku.

16


Rastriginova funkce, multimodáln´ı, separabiln´ı, d ∑ [ 2 ] f (x) = 10d + xi − 10 cos(2πxi ) i=1

xi ∈ [−5.12, 5.12], x∗ = (0, 0, . . . 0), f (x∗ ) = 0. Nalezen´ı minima této funkce je povaˇzováno za obt´ıˇznou u ´lohu globáln´ı optimalizace. Graf této funkce pro d = 2 je na následuj´ıc´ım obrázku.


17

Schwefelova funkce, multimodáln´ı, separabiln´ı, je zde uvedena ve tvaru s aditivn´ım ˇclenem za rovn´ıtkem, kter´ y posouvá hodnoty funkce tak, aby jej´ı hodnota v globáln´ım minimu byla rovna nule stejnˇe jako u ostatn´ıch zde uveden´ ych testovac´ıch funkc´ı. d ) (√ ∑ | xi | , f (x) = 418.9828872724338 d − xi sin i=1

. . xi ∈ [−500, 500], x∗ = (s, s, . . . s), s = 420.968746, f (x∗ ) = 0, pro d < 80 plat´ı, ˇze 0 <= f (x∗ ) < 1e − 10. Nalezen´ı minima této funkce je povaˇzováno za stˇrednˇe obt´ıˇznou u ´lohu globáln´ı optimalizace, jej´ı globáln´ı minimum je v prohledávaném prostoru velmi vzdálené od druhého nejlepˇs´ıho lokáln´ıho minima. Graf této funkce pro d = 2 (bez aditivn´ıho ˇclenu posouvaj´ıc´ıho funkˇcn´ı hodnoty) je na následuj´ıc´ım obrázku.

Nˇekteˇr´ı autoˇri oprávnˇenˇe povaˇzuj´ı funkce, u kter´ ych je globáln´ı minimum ve stˇredu prohledávané domény (z uveden´ ych funkc´ı je to prvn´ı De Jongova, Ackleyho, Griewankova a Rastriginova funkce), tj. v tˇeˇziˇsti d-rozmˇerného kvádru D, za nevhodné pro testován´ı stochastick´ ych algoritm˚ u, nebot’ u nich lze naj´ıt globáln´ı minimum pouh´ ym pr˚ umˇerován´ım náhodnˇe vygenerovan´ ych bod˚ u v D. Tento nedostatek testovac´ıch funkc´ı s bodem globáln´ıho minima v tˇeˇziˇsti domény D lze snadno odstranit. M´ısto vyhodnocen´ı funkce v bodˇe x se vyhodnocuje funkce v bodˇe x − o. Bod o ∈ D znamená posun (shift). Globáln´ı minimum takové posunuté funkce je pak x∗ = o. Posun o m˚ uˇzeme zadat pˇredem nebo generovat náhodnˇe pˇred kaˇzd´ ym spuˇstˇen´ım optimalizace testovac´ı funkce. Implementace tˇechto testovac´ıch funkc´ı je v Matlabu velmi jednoduchá, jak vid´ıte v následuj´ıc´ıch zdrojov´ ych textech.

18


function y=dejong1(x) % First Dejong y=sum(x.*x); function y=rosen(x) % Rosenbrockova funkce n=length(x); a=(x(1:n-1).^2-x(2:n)).^2; b=1-x(1:n-1); y=sum(100*a+b.^2); function y=ackley(x) % Ackley’s function, <-30,30>, glob. min f(0)=0 d=length(x); a=-20*exp(-0.02*sqrt(sum(x.*x))); b=-exp(sum(cos(2*pi*ones(1,d).*x))/d); y=a+b+20+exp(1); function y=griewank(x) % Griewank’s function, <-400,400>, glob. min f(0)=0 d=length(x); a=sum(x.*x)/4000; j=1:d; b=prod(cos(x./sqrt(j))); y=a-b+1; % function y=rastrig(x) % Rastrigin’s function, <-5.12,5.12>, glob. min f(0)=0 d=length(x); sz=size(x); if sz(1)==1 x=x’; end y=10*d+sum(x.*x-10*cos(2*pi*x)); %


19

function y=schwefel0(x) % Schwefel’s function, <-500,500>, % glob. min f(x_star)~0, (for d<80 0<=f(x_star)<1e-10), % x_star=[s,s,...,s], s=420.968746 d=length(x); sz=size(x); if sz(1)==1 x=x’; end y = 418.9828872724338*d - sum(x.*sin(sqrt(abs(x)))); %

Shrnut´ı: - Testovac´ı funkce, numerické testován´ı algoritm˚ u - R˚ uzné typy ukonˇcen´ı prohledáván´ı ˇ - Casov´ a nároˇcnost mˇeˇrená poˇctem vyhodnocen´ı u ´ˇcelové funkce, spolehlivost nalezen´ı globáln´ıho minima - Vyhodnocen´ı experimentu

Kontroln´ı ot´ azky: 1. Proˇc je nutné posuzovat algoritmy nejen z jednoho bˇehu, ale z v´ıce opakován´ı? 2. Jak popsat specifikaci testovac´ıho problému? 3. Jak specifikovat testovan´ y algoritmus?

´ NAHODN ´ ´ PROHLEDAV ´ AN ´ Í 4 SLEPE E

20

4

Slep´ e n´ ahodn´ e prohled´ av´ an´ı

Pr˚ uvodce studiem: V této kapitole se seznám´ıte s nejjednoduˇsˇs´ım stochastick´ ym algoritmem pro hledán´ı globáln´ıho minima. Kromˇe toho uvid´ıte i snadnost zápisu takov´ ych algoritm˚ u v Matlabu a seznám´ıte se s pojmem konvergence stochastického algoritmu. Poˇc´ıtejte asi se dvˇema hodinami studia.

Náhodné prohledáván´ı se také nˇekdy naz´ yvá slep´ y algoritmus, nebot’ se v nˇem opa∏ kovanˇe generuje náhodné ˇreˇsen´ı (nov´ y bod y v souvislé oblasti D = di=1 ⟨ai , bi ⟩ z rovnomˇerného spojitého rozdˇelen´ı) a zapamatovává se tehdy, kdyˇz je lepˇs´ı neˇz pˇredt´ım nalezené ˇreˇsen´ı. Náhodné prohledáván´ı je pˇr´ıklad velmi jednoduchého stochastického algoritmu pro hledán´ı globáln´ıho minima vyuˇz´ıvaj´ıc´ıho velmi prostou heuristiku, kterou bychom mohli slovnˇe vyjádˇrit jako “zkus bez rozmyslu cokoliv, co nen´ı zakázáno”. Slep´ y algoritmus lze zapsat v Matlabu na pár ˇrádc´ıch: function [x,fx] = blindsearch(fn_name, a, b, t) % input parameters: % fn_name function to be minized (M file) % a, b row vectors, limits of search space % t iteration number % output: % fx minimal function value found by the search % x minimum point found by the search d = length(a); fx = realmax; for i=1:t y = a + (b - a).*rand(1,d); fy = feval(fn_name, y); if fy < fx x = y; fx = fy; end end

21

U tohoto algoritmu lze také pomˇernˇe jednoduˇse dokázat jeho teoretickou konvergenci. Plat´ı, ˇze s rostouc´ım poˇctem iterac´ı se nalezené ˇreˇsen´ı x pˇribliˇzuje globáln´ımu minimu x∗ : lim P (∥ (x − x∗ ) ∥< ε | t) = 1 ε > 0, (4.1) t→∞

∗

kde ∥ (x − x ) ∥ je vzdálenost nalezeného ˇreˇsen´ı x od globáln´ıho minima x∗ (norma vektoru x − x∗ ) a P (∥ (x − x∗ ) ∥< ε | t) je pravdˇepodobnost jevu (∥ (x − x∗ ) ∥< ε) za podm´ınky, ˇze bylo provedeno t iterac´ı. K tvrzen´ı vyjádˇrenému rovnic´ı (4.1) dojdeme následuj´ıc´ı u ´vahou: Vˇsechny body ∗ splˇ nuj´ıc´ı podm´ınku ∥ (x − x ) ∥< ε jsou vnitˇrn´ımi body d-rozmˇerné koule Dε o polomˇeru ε). Jelikoˇz ε > 0, je“objem”této koule kladn´ y, pˇresnˇeji m´ıra této mnoˇziny λ(Dε ) > 0. Oznaˇcme dále m´ıru mnoˇziny D jako λ(D). Generujeme-li bod y ∈ D z rovnomˇerného rozdˇelen´ı, pak pravdˇepodobnost jevu A ≡ {y ∈ Dε } je P (A) =

λ(Dε ) λ(D)

a pravdˇepodobnost jevu opaˇcného je ¯ = 1 − λ(Dε ) P (A) λ(D) Pravdˇepodobnost, ˇze jev A¯ nastane v t po sobˇe jdouc´ıch opakován´ıch je ( )t λ(D ) ε P (A¯ | t) = 1 − . λ(D) ( Jelikoˇz 0 < 1 −

λ(Dε ) λ(D)

) < 1, je potom lim P (A¯ | t) =

t→∞

( )t λ(Dε ) 1− = 0, λ(D)

pravdˇepodobnost jevu opaˇcného je rovna 1, a tedy rovnice (4.1) plat´ı. Lze ukázat, ˇze nejenom slepé prohledáván´ı, ale dokonce kaˇzd´ y stochastick´ y algoritmus, ve kterém se jin´ y zp˚ usob generován´ı nového bodu y stˇr´ıdá s jeho slep´ ym generován´ım z rovnomˇerného rozdˇelen´ı v prohledávaném prostoru D a toto slepé generován´ı má v kaˇzdém kroku kladnou pravdˇepodobnost (pˇresnˇeji, kdyˇz tyto pravdˇepodobnosti pro jednotlivé kroky algoritmu tvoˇr´ı divergentn´ı ˇradu, ˇcehoˇz lze snadno dosáhnout volbou konstantn´ı pravdˇepodobnosti v kaˇzdém kroku), tak takov´ y kombinovan´ y algoritmus je konvergentn´ı ve smyslu rov. (4.1). Bohuˇzel rovnice (4.1) nám neposkytuje ˇza´dnou informaci o kvalitˇe nalezeného ˇreˇsen´ı po koneˇcném poˇctu iterac´ı. Jak uvid´ıme pozdˇeji, s podobn´ ym problémem se budeme setkávat i u dalˇs´ıch stochastick´ ych algoritm˚ u. Teoretická konvergence podle

´ NAHODN ´ ´ PROHLEDAV ´ AN ´ Í 4 SLEPE E

22

rovnice (4.1) je z praktického hlediska velmi slabé tvrzen´ı a pro posouzen´ı konvergence stochastick´ ych algoritm˚ u a porovnán´ı jejich efektivity jsou d˚ uleˇzitˇejˇs´ı experimentáln´ı v´ ysledky na testovac´ıch funkc´ıch neˇz teoretick´ y d˚ ukaz konvergence ve formˇe rovnice (4.1). Shrnut´ı: - Generován´ı bodu z rovnomˇerného rozdˇelen´ı na D - Konvergence stochastického algoritmu

Kontroln´ı ot´ azky: 1. Jakou heuristiku uˇz´ıvá slep´ y algoritmus pro generován´ı nového bodu? 2. Je ve slepém náhodném prohledáván´ı modelován nˇejak´ y proces uˇcen´ı?

23

5

ˇ ızen´ R´ e n´ ahodn´ e prohled´ av´ an´ı

Pr˚ uvodce studiem: V této kapitole se poprvé seznám´ıte s algoritmem, kter´ y pracuje s populac´ı jedinc˚ u, tj. bod˚ u (kandidát˚ u ˇreˇsen´ı) v D. Algoritmus ˇr´ızeného náhodného prohledáván´ı je jednoduch´ y a efektivn´ı stochastick´ y algoritmus pro hledán´ı globáln´ıho minima. Pochopen´ı tohoto algoritmu je dobr´ ym odrazov´ ym m˚ ustkem k evoluˇcn´ım algoritm˚ um. Poˇc´ıtejte asi se dvˇema hodinami studia a nˇekolika hodinami vˇenovan´ ymi implementaci algoritmu v Matlabu a jeho ovˇeˇren´ı na testovac´ıch funkc´ıch. ˇ ızené náhodné prohledáván´ı (controlled random search, CRS) je pˇr´ıkladem velmi R´ jednoduchého a pˇritom efektivn´ıho stochastického algoritmu (ˇcili heuristiky) pro hledán´ı minima v souvislé oblasti D. Prvn´ı verzi algoritmu CRS navrhl Price [18] v sedmdesát´ ych letech minulého stolet´ı. Pracuje se s populac´ı N bod˚ u – kandidát˚ u ˇreˇsen´ı v prohledávaném prostoru D a z nich se pomoc´ı nˇejaké lok´ aln´ı heuristiky generuje nov´ y bod y, kter´ y m˚ uˇze b´ yt zaˇrazen do populace m´ısto dosud nejhorˇs´ıho bodu. Poˇcet vygenerovan´ ych bod˚ u populace N je vˇetˇs´ı neˇz dimenze d prohledávaného prostoru D. Price jako lokáln´ı heuristiku pro generován´ı nového bodu y uˇz´ıval reflexi simplexu, která je známa z velmi populárn´ı simplexové metody hledán´ı minima, kterou navrhli pˇred zhruba 50 léty Nelder a Mead [16]. Simplexová metoda je vyuˇzita také ve vestavˇené funkci fminsearch v Matlabu jako prostˇredek pro hledán´ı minima zadané funkce. Reflexe simplexu v CRS se realizuje takto: Z populace o velikosti N , N > d se náhodnˇe vybere d + 1 bod˚ u tvoˇr´ıc´ıch simplex (optimisticky spoléháme na to, ˇze tyto náhodnˇe vybrané body jsou nekomplanárn´ı, tj. libovoln´ ych d vektor˚ u z tˇechto d+1 vektor˚ u tvoˇr´ı bázi d-rozmˇerného prostoru). Pak bod simplexu s nejvˇetˇs´ı funkˇcn´ı hodnotou pˇreklop´ıme kolem tˇeˇziˇstˇe zb´ yvaj´ıc´ıch bod˚ u simplexu a z´ıskáme tak nového potenciáln´ıho kandidáta ˇreˇsen´ı y. Reflexe v simplexu je vyjádˇrena vztahem y = g + (g − xH ) = 2g − xH ,

(5.1)

kde xH je bod simplexu s nejvˇetˇs´ı hodnotou u ´ˇcelové funkce a g je tˇeˇziˇstˇe zb´ yvaj´ıc´ıch d bod˚ u simplexu, jehoˇz souˇradnice spoˇc´ıtáme jako pr˚ umˇery tˇechto d bod˚ u simplexu. Graficky je reflexe znázornˇena na následuj´ıc´ım obrázku.

ˇ ÍZENE ´ NAHODN ´ ´ PROHLEDAV ´ AN ´ Í 5 R E

24

8 7

x

L

6

x

H

5

g

4

y 3 2 1 0

0

1

2

3

4

5

6

7

8

Pokud v novém bodu y je funkˇcn´ı hodnota f (y) menˇs´ı neˇz v nejhorˇs´ım bodu celé populace, pak je tento nejhorˇs´ı bod nahrazen bodem y. Nahrazen´ım nejhorˇs´ıho bodu populace nov´ ym bodem y dosahujeme toho, ˇze populace se koncentruje v okol´ı dosud nalezeného bodu s nejmenˇs´ı funkˇcn´ı hodnotou. Algoritmus m˚ uˇzeme jednoduˇse zapsat ve strukturovaném pseudokódu následuj´ıc´ım zp˚ usobem: generuj populaci P , tj. N bod˚ u náhodnˇe v D repeat najdi xworst ∈ P takové, ˇze f (xworst ) ≥ f (x), repeat vyber z P simplex y := reflexe simplexu, y ∈ D until f (y) < f (xworst ); xworst := y; until podm´ınka ukonˇcen´ı;

x∈P

Reflexe simplexu podle rovnice (5.1) samozˇrejmˇe nen´ı jediná moˇznost, jak v algoritmu CRS generovat nov´ y bod. Price pozdˇeji navrhl jinou lokáln´ı heuristiku, kdy nov´ y bod y se generuje podle (5.1), ale do simplexu je vˇzdy zaˇrazen nejlepˇs´ı bod populace xmin s funkˇcn´ı hodnotou fmin , fmin ≤ f (x), x ∈ P a zb´ yvaj´ıc´ıch d bod˚ u simplexu se pak vybere náhodnˇe z ostatn´ıch bod˚ u populace. Dalˇs´ı moˇznost je znáhodnˇená reflexe, která byla navrˇzena v 90. létech. Znáhodnˇená reflexe v simplexu je popsána vztahem y = g + U (g − xH ),

(5.2)

25

U je náhodná veliˇcina vhodného rozdˇelen´ı. Graficky je taková reflexe znázornˇena na následuj´ıc´ım obrázku.

10 9 D 8 7 6 U(g−x) 5

x

g

2

4

H

y

4 3 2 1 0 −1

0

6

8

10

Pokud bychom uˇzili reflexi podle (5.1), byla by vzdálenost bod˚ u g a y byla stejná jako vzdálenost bod˚ u x a g, tj. hodnota U v rovnici (5.2) by byla konstantn´ı bˇehem prohledáván´ı a rovna jedné. Pˇri uˇzit´ı znáhodnˇené reflexe je tato vzdálenost urˇcována náhodnˇe podle zvoleného rozdˇelen´ı náhodné veliˇciny U . Na testovac´ıch u ´lohách bylo ovˇeˇrováno nˇekolik rozdˇelen´ı náhodné veliˇciny U . Osvˇedˇcilo se rovnomˇerné rozdˇelen´ı na [0, α), kde α je vstupn´ı parametr algoritmu. Podle empirick´ ych v´ ysledk˚ u test˚ u nejrychleji algoritmus konvergoval na vˇetˇsinˇe testovac´ıch funkc´ı pˇri hodnotách α ≈ 4. Stˇredn´ı hodnota je E U = α/2. Za pozornost stoj´ı uˇzit´ı rovnomˇerného rozdˇelen´ı náhodné veliˇciny U na [s, α − s), kde α > 0 a s ∈ (0, α/2) jsou vstupn´ı parametry této lokáln´ı heuristiky. Stˇredn´ı hodnota je opˇet E U = α/2. Je vˇsak moˇzno uˇz´ıt i jiná rozdˇelen´ı produkuj´ıc´ı nezáporné hodnoty náhodné veliˇciny U , napˇr. lognormáln´ı rozdˇelen´ı, pˇr´ıpadnˇe i jiné lokáln´ı heuristiky pro generován´ı nového bodu y, neˇz je reflexe. Reflexe podle rovnic (5.1) nebo (5.2) nezaruˇcuje, ˇze novˇe vygenerovan´ y bod y bude v prohledávaném prostoru D. Pokud nastane situace, ˇze y ̸∈ D, jsou r˚ uzné moˇznosti, jak postupovat dále. Bud’ generujeme nov´ y bod opakovanˇe tak dlouho, neˇz se po∏d ypoˇcetnˇe u ´spornˇejˇs´ı postup. daˇr´ı vygenerovat y ∈ D = i=1 ⟨ai , bi ⟩ nebo pouˇzijeme v´ Napˇr. je moˇzné vygenerovat souˇradnice yi ̸∈ ⟨ai , bi ⟩ náhodnˇe nebo uˇz´ıt tzv. zrcadlen´ı (perturbaci), kdy ty souˇradnice yi ̸∈ ⟨ai , bi ⟩ pˇreklop´ıme dovnitˇr prohledávaného prostoru D kolem pˇr´ısluˇsné strany d-rozmˇerného kvádru D. Toto zrcadlen´ı znázorˇ nuje následuj´ıc´ı obrázek, jeho algoritmus je zapsán jako funkce zrcad v Matlabu.


26

10 9 D 8 7 6 y mimo D

y po perturbaci

5 4 3 2 1 0 −1 −1

0

1

2

3

4

5

6

7

8

9

10

% zrcadleni, perturbation y into function y = zrcad(y, a, b) zrc = find(y < a | y > b); % najdi pretekajici dimenze for i = zrc % preklop while (y(i) < a(i) || y(i) > b(i)) if y(i) > b(i) y(i) = 2*b(i) - y(i); elseif y(i) < a(i) y(i) = 2*a(i) - y(i); end end end Pˇri stochastickém prohledáván´ı podm´ınku ukonˇcen´ı formuluje uˇzivatel. Obvykle je moˇzné uˇz´ıt podm´ınku ukonˇcen´ı formulovanou tak, ˇze nejlepˇs´ı a nejhorˇs´ı bod populace se ve funkˇcn´ıch hodnotách liˇs´ı jen málo, tedy f (xworst ) − f (xmin ) < ε,

(5.3)

kde ε > 0 je vstupn´ı parametr algoritmu. Jinou moˇznost´ı je hledán´ı ukonˇcit pˇri dosaˇzen´ı zadaného poˇctu vyhodnocen´ı funkce, pˇr´ıpadnˇe tyto podm´ınky kombinovat – viz kap. 3.

27

Algoritmus CRS m˚ uˇzeme celkem rychle napsat v Matlabu takto: function [x_star, fn_star, func_evals]=... crs1(fn_name, a, b, N, my_eps, max_evals, alfa, shift) d = length(a); P = zeros(N,d+1); for i = 1:N P(i,1:d) = a + (b - a).*rand(1,d); P(i,d+1) = feval(fn_name,P(i,1:d) - shift); end % 0-th generation initialized [fmax, indmax] = max(P(:,d+1)); [fmin, indmin] = min(P(:,d+1)); func_evals = N; while (fmax - fmin > my_eps) && (func_evals < d*max_evals)% main loop y = reflrand(P, alfa); y = zrcad(y, a, b); % perturbation fy = feval(fn_name,y - shift); func_evals = func_evals + 1; if fy < fmax % y is good P(indmax,:) = [y, fy]; [fmax, indmax] = max(P(:,d+1)); [fmin, indmin] = min(P(:,d+1)); end end % main loop - end x_star = P(indmin,1:d); fn_star = fmin; end Funkce crs1 kromˇe zrcadlen´ı zrcad volá dalˇs´ı dvˇe funkce, a to reflrand, která realizuje znáhodnˇenou reflexi podle vztahu (5.2) s U rovnomˇernˇe rozdˇeleném na [0, α) a funkci nahvyb_expt, která vyb´ırá z N bod˚ u populace náhodnˇe d + 1 bod˚ u do simplexu. Prostudujte pozornˇe implementaci vˇsech tˇechto funkc´ı, nebot’ se v nich m˚ uˇzete seznámit i s r˚ uzn´ ymi technikami efektivn´ıho vyuˇzit´ı moˇznost´ı Matlabu, které znaˇcnˇe sniˇzuj´ı délku zdrojového kódu a t´ım i pracnost implementace. Funkce nahvyb_expt a zrcad vyuˇzijeme i pozdˇeji v nˇekter´ ych evoluˇcn´ıch algoritmech.

28


function y = reflrand(P, alpha) % znahodnena reflexe simplexu N = length(P(:,1)); d = length(P(1,:)) - 1; vyb = nahvyb_expt(N, d+1); S = P(vyb,:); % simplex S vybran nahodne z P [x, indx] = max(S(:, d+1)); % nejhorsi bod v S x = S(indx, 1:d); S(indx,:)= []; % zrus radek s nejhorsi hodnotou S(:,d+1) = []; % zrus posledni sloupec s funkcni hodnotou g = mean(S); % teziste y = g + alpha * rand(1) * (g - x); % random sample, k of N without repetition, % numbers given in vector expt are not included function vyb = nahvyb_expt(N, k, expt) opora = 1:N; if nargin == 3 opora(expt) = []; end vyb = zeros(1,k); for i = 1:k index = 1 + fix(rand(1)*length(opora)); vyb(i) = opora(index); opora(index) = []; end Algoritmus CRS je pro nás prvn´ım algoritmem, kter´ y pracuje s populac´ı bod˚ u (potenciáln´ıch kandidát˚ u moˇzného hledaného ˇreˇsen´ı) a tuto populaci nechává vyv´ıjet v pr˚ ubˇehu prohledávac´ıho procesu. Proto nˇekteˇr´ı autoˇri algoritmus CRS ˇrad´ı mezi evoluˇcn´ı algoritmy, i kdyˇz z princip˚ u evoluˇcn´ıch algoritm˚ u jsou uplatnˇeny pouze nˇekteré, jak uvid´ıme v následuj´ıc´ıch kapitolách.

29

Shrnut´ı: - Základn´ı principy algoritmu CRS - Populace bod˚ u náhodnˇe vygenerovan´ ych v D - Lokáln´ı heuristika generuj´ıc´ı nov´ y pokusn´ y bod - Simplex náhodnˇe vybran´ y z populace - Reflexe a znáhodnˇená reflexe - Podm´ınka ukonˇcen´ı

Kontroln´ı ot´ azky: 1. Jak urˇc´ıte souˇradnice tˇeˇziˇstˇe n bod˚ u v D? 2. Co je to zrcadlen´ı a kdy se uˇz´ıvá? 3. Jak lze zformulovat podm´ınku ukonˇcen´ı? 4. Jak byste experimentálnˇe porovnali u ´ˇcinnost algoritmu CRS a slepého prohledáván´ı na vybran´ ych testovac´ıch u ´lohách?

ˇ Í ALGORITMY 6 EVOLUCN

30

6

Evoluˇ cn´ı algoritmy

Pr˚ uvodce studiem: V této kapitole se struˇcnˇe seznám´ıme s histori´ı a principy evoluˇcn´ıch algoritm˚ u. Na tuto krátkou kapitolu poˇc´ıtejte asi s hodinou studia.

V posledn´ıch desetilet´ıch se s pomˇern´ ym u ´spˇechem pro hledán´ı globáln´ıho minima funkc´ı uˇz´ıvaj´ı stochastické algoritmy zejména evoluˇcn´ıho typu. Podrobn´ y popis této problematiky naleznete v knihách Goldgerga [8], Michalewicze [13] nebo Bäcka [2]. Rozvoj evoluˇcn´ıch algoritm˚ u je záleˇzitost´ı posledn´ıch desetilet´ı a je podm´ınˇen rozvojem poˇc´ıtaˇc˚ u a pokroky v informatice. Pˇrelomov´ ymi pracemi zavádˇej´ıc´ımi biologickou terminologii do model˚ u hledán´ı globáln´ıho extrému jsou ˇclánky Schwefela a Rechenberga ze ˇsedesát´ ych let minulého stolet´ı o evoluˇcn´ı strategii, práce L. Fogela o evoluˇcn´ım programován´ı a Hollandova kniha o genetick´ ych algoritmech [9]. T´ım byl odstartován bouˇrliv´ y rozsáhl´ y rozvoj evoluˇcn´ıch algoritm˚ u a jejich aplikac´ı. Evoluˇcn´ı algoritmy jsou ve své podstatˇe jednoduch´ ymi modely Darwinovy evoluˇcn´ı teorie v´ yvoje populac´ı. Charakteristické pro nˇe je to, ˇze pracuj´ı s populac´ı (tvoˇrenou jedinci – kandidáty moˇzn´ ych ˇreˇsen´ı), která se v pr˚ ubˇehu hledán´ı vyv´ıj´ı (vznikaj´ı nové generace). Tento v´ yvoj se uskuteˇcn ˇuje aplikac´ı evoluˇcn´ıch oper´ ator˚ u, nejd˚ uleˇzitˇejˇs´ı evoluˇcn´ı operátory jsou tyto: • selekce – silnˇejˇs´ı jedinci z populace maj´ı vˇetˇs´ı pravdˇepodobnost svého pˇreˇzit´ı i pˇredán´ı sv´ ych vlastnost´ı potomk˚ um, • kˇr´ıˇzen´ı (rekombinace) – dva nebo v´ıce jedinc˚ u z populace si vymˇen´ı informace a vzniknou tak nov´ı jedinci kombinuj´ıc´ı vlastnosti rodiˇc˚ u, • mutace – genetická informace zakódovaná v jedinci m˚ uˇze b´ yt náhodnˇe modifikována. Dále se mezi evoluˇcn´ı operátory poˇc´ıtá migrace, která se uˇz´ıvá v tzv. paraleln´ıch evoluˇcn´ıch algoritmech, kdy se modeluje v´ yvoj v´ıce populac´ı vedle sebe a jejich vzájemné ovlivˇ nován´ı, ale takov´ ymi algoritmy se v tomto kurzu nebudeme zab´ yvat. Evoluˇcn´ı algoritmy jsou heuristické postupy, které nˇejak´ ym zp˚ usobem modifikuj´ı populaci tak, aby se jej´ı vlastnosti zlepˇsovaly. O nˇekter´ ych tˇr´ıdách evoluˇcn´ıch algoritm˚ u je dokázáno, ˇze nejlepˇs´ı jedinci populace se skuteˇcnˇe pˇribliˇzuj´ı ke globáln´ımu extrému. Evoluˇcn´ı algoritmy byly a jsou pˇredmˇetem intenzivn´ıho v´ yzkumu a poˇcet publikac´ı z této oblasti je velk´ y. Jedn´ım z hlavn´ıch motiv˚ u jsou pˇredevˇs´ım aplikace v praktick´ ych problémech, které jin´ ymi metodami nejsou ˇreˇsitelné. Dalˇs´ımi motivy pro roz-

31

voj evoluˇcn´ıch algoritm˚ u jsou v´ yzkum umˇelé inteligence a teorie uˇcen´ı. Tento rozvoj nezastavil ani tzv. ”No Free Lunch Theorem”[27], kter´ y ukazuje, ˇze nelze naj´ıt univerzálnˇe nejlepˇs´ı stochastick´ y algoritmus pro globáln´ı optimalizaci. Porovnáváme-li dva algoritmy, vˇzdy lze naj´ıt skupinu problém˚ u, ve kter´ ych bude prvn´ı algoritmus u ´ˇcinnˇejˇs´ı a jinou skupinu problém˚ u, ve kter´ ych bude u ´ˇcinnˇejˇs´ı algoritmus druh´ y.

Shrnut´ı: - Populace a jej´ı v´ yvoj - Evoluˇcn´ı operátory - No Free Lunch Theorem

Kontroln´ı ot´ azky: ˇ ım jsou evoluˇcn´ı algoritmy inspirovány? 1. C´ 2. Proˇc byly evoluˇcn´ı algoritmy navrˇzeny a proˇc jsou dále zkoumány a rozv´ıjeny?


32

6.1

Genetick´ e algoritmy

Pr˚ uvodce studiem: Tato struˇcná kapitola má poslouˇzit k základn´ı orientaci v principech genetick´ ych algoritm˚ u a upozornit na postupy, které jsou inspirativn´ı pro jiné evoluˇcn´ı algoritmy. Poˇc´ıtejte asi se dvˇema hodinami studia.

Genetické algoritmy od svého vzniku v 70. letech [9] pˇredstavuj´ı mohutn´ y zdroj inspirace pro dalˇs´ı evoluˇcn´ı algoritmy, pro zkoumán´ı umˇelé inteligence a pro rozvoj soft computingu. O genetick´ ych algoritmech existuje ˇrada specializovan´ ych monografi´ı, viz napˇr. Goldberger [8], Michalewicz [13] nebo Bäck [2]. Velmi dobˇre a srozumitelnˇe napsan´ yu ´vod do genetick´ ych algoritm˚ u najdete v dostupné knize Kvasniˇcky, Posp´ıchala a Tiˇ na [12]. Zde se omez´ıme jen na struˇcné vysvˇetlen´ı základn´ıch princip˚ u tˇechto algoritm˚ u. Základn´ı myˇslenkou genetick´ ych algoritm˚ u je analogie s evoluˇcn´ımi procesy prob´ıhaj´ıc´ımi v biologick´ ych systémech. Podle Darwinovy teorie pˇrirozeného v´ ybˇeru pˇreˇz´ıvaj´ı jen nejlépe pˇrizp˚ usoben´ı jedinci populace. M´ırou pˇrizp˚ usoben´ı je tzv. “fitness” jedince. V biologii je fitness chápána jako relativn´ı schopnost pˇreˇzit´ı a reprodukce genotypu jedince. Biologická evoluce je zmˇena obsahu genetické informace populace v pr˚ ubˇehu mnoha generac´ı smˇerem k vyˇsˇs´ım hodnotám fitness. Jedinci s vyˇsˇs´ı fitness maj´ı vˇetˇs´ı pravdˇepodobnost pˇreˇzit´ı a vˇetˇs´ı pravdˇepodobnost reprodukce sv´ ych gen˚ u do generace potomk˚ u. Kromˇe reprodukce se v populaˇcn´ım v´ yvoji uplatˇ nuje i tzv. mutace, coˇz je náhodná zmˇena genetické informace nˇekter´ ych jedinc˚ u v populaci. V genetick´ ych algoritmech je fitness kladné ˇc´ıslo pˇriˇrazené genetické informaci jedince. Tato genetická informace jedince (chromozóm) se obvykle vyjadˇruje bitov´ ym ˇretˇezcem. Populaci jedinc˚ u pak modelujeme jako populaci chromozóm˚ u a kaˇzdému chromozómu (bitovému ˇretˇezci) um´ıme pˇriˇradit hodnotu u ´ˇcelové funkce a podle vhodného pˇredpisu um´ıme pˇriˇradit i jeho fitness, podrobnˇeji viz [12]. Obvykle je chromozómem binárn´ı vektor konstantn´ı délky k α = (α1 , α2 , . . . , αk ) ∈ {0, 1}k a populace velikosti N je potom mnoˇzina takov´ ych chromozóm˚ u P = {α1 , α2 , . . . , αN }

6.1 Genetické algoritmy

33

Pokud hledáme globáln´ı minimum u ´ˇcelové funkce f , pak fitness F je nˇejaké zobrazen´ı, které vyhovuje následuj´ıc´ı podm´ınce f (αi ) ≤ f (αj ) ⇒ F (αi ) ≥ F (αj ) > 0, i, j = 1, 2, . . . , N Nejjednoduˇsˇs´ı zp˚ usob, jak vyhovˇet této podm´ınce, je lineárn´ı zobrazen´ı funkˇcn´ıch hodnot na fitness, kdy hodnotu fitness urˇc´ıme podle vztahu F (α) =

Fmax − Fmin fmin Fmin − fmax Fmax f (α) + , fmin − fmax fmin − fmax

kde fmin a fmax je nejmenˇs´ı a nejvˇetˇs´ı hodnota funkce v populaci, Fmin a Fmax je minimáln´ı a maximáln´ı hodnota fitness, které obvykle vol´ıme Fmax = 1 a Fmin = ϵ, kde ϵ je malé kladné ˇc´ıslo, napˇr. ϵ = 0.01. Pak vztah pro fitness má tvar F (α) =

1 [(1 − ϵ)f (α) + ϵfmin − fmax ] . fmin − fmax

Hodnoty fitness je vhodné renormalizovat, tj. pˇrepoˇc´ıtat tak, aby jejich souˇcet byl roven jedné. Renormalizovaná fitness i-tého jedince je F (αi ) F ′ (αi ) = ∑N j=1 F (αj ) a tato hodnota pak udává i pravdˇepodobnost pˇreˇzit´ı jedince a jeho u ´ˇcasti na reprodukci své genetické informace, tj. vytvoˇren´ı potomka. V´ ybˇer (selekci) jedince s prav′ dˇepodobnost´ı F (αi ) ukazuje následuj´ıc´ı algoritmus naz´ yvan´ y také ruleta. Vstupn´ım parametrem fits je vektor hodnot fitness vˇsech jedinc˚ u populace, hodnoty fitness nemus´ı b´ yt renormalizovány. Funkce vrac´ı ˇc´ıslo (index) vybraného jedince, pravdˇepodobnost v´ ybˇeru je u ´mˇerná hodnotˇe fitness.


34

function res=roulette_simple(fits) % % returns an integer from [1, length(fits)] % with probability proportional % to fits(i)/ sum fits h = length(fits); ss = sum(fits); cp = zeros(1, h); cp(1) = fits(1); for i=2:h cp(i) = cp(i-1) + fits(i); end cp = cp/ss; res = 1 + fix(sum(cp < rand(1)));

Reprodukce prob´ıhá kˇr´ıˇzen´ım dvou chromozóm˚ u tak, ˇze si vymˇen´ı své geny (ˇcásti bitov´ ych ˇretˇezc˚ u). Ze dvou rodiˇc˚ u α, β v jednobodovém kˇr´ıˇzen´ı vzniknou dva potomci podle následuj´ıc´ıho schématu: (α1 , . . . , αl , α(l+1) , . . . , αk ) ↘ (α1 , . . . , αl , β(l+1) , . . . , βk ) (β1 , . . . , βl , β(l+1) , . . . , βk ) ↗ (β1 , . . . , βl , α(l+1) , . . . , αk ) Bod kˇr´ıˇzen´ı l se vol´ı náhodnˇe. Nˇekdy se také uˇz´ıvá dvoubodové kˇr´ıˇzen´ı, ve kterém se urˇc´ı náhodnˇe dvˇe pozice pro kˇr´ıˇzen´ı, l1 , l2 , l2 > l1 a v bitov´ ych ˇretˇezc´ıch se vymˇen´ı u ´seky αl1 , . . . , αl2 ↔ βl1 , . . . , βl2 Mutace je vˇetˇsinou v genetick´ ych algoritmech implementována jako zmˇena hodnoty náhodnˇe vybraného bitu v chromozómu, tj. hodnota 0 se zmˇen´ı na 1, hodnota 1 na 0. Mutac´ı je zajiˇstˇeno, aby mohla v populaci vzniknout genetická informace, která v n´ı v pˇredcházej´ıc´ıch generac´ı nebyla, nebo se obnovila genetické informace ztracené v pr˚ ubˇehu dosavadn´ıho v´ yvoje. V d˚ usledku toho pak algoritmus m˚ uˇze uniknout z oblasti lokáln´ıho minima, ke kterému by smˇeˇroval, pokud bychom uˇz´ıvali pouze kˇr´ıˇzen´ı. Naopak, pokud bychom uˇz´ıvali pouze mutaci, genetick´ y algoritmus by se choval podobnˇe jako slepé náhodné prohledáván´ı. Lze ukázat, ˇze genetick´ y algoritmus konverguje ke globáln´ımu extrému. Genetické algoritmy na rozd´ıl od ostatn´ıch algoritm˚ u uvádˇen´ ych v tomto textu reprezentuj´ı jedince jako bitov´ y ˇretˇezec, nikoliv jako vektor hodnot v pohyblivé ˇrádové ˇca´rce. To je ˇcasto v´ yhodou pˇri ˇreˇsen´ı tzv. diskrétn´ıch u ´loh globáln´ı optimalizace (kombinatorické u ´lohy jako je problém obchodn´ıho cestuj´ıc´ıho), ale pˇrináˇs´ı to kom-

6.1 Genetické algoritmy

35

plikace v problémech, kdy je prohledávaná oblast D moˇzn´ ych ˇreˇsen´ı spojitá. Pokud bychom u ´seky chromozómu interpretovali jako celá ˇc´ısla (datov´ y typ integer), pak se pot´ ykáme s obt´ıˇz´ı, ˇze sousedn´ı ˇc´ıselné hodnoty jsou reprezentovány ˇretˇezci, které se mohou podstatnˇe liˇsit, dokonce i ve vˇsech bitech, napˇr. dekadická hodnota 7 je binárnˇe (0111), následuj´ıc´ı ˇc´ıselná hodnota 8 je (1000). Tuto nepˇr´ıjemnost je nutno v genetick´ ych algoritmech nˇejak oˇsetˇrit, pokud maj´ı b´ yt u ´spˇeˇsnˇe pouˇz´ıvány i pro ˇreˇsen´ı spojit´ ych problém˚ u globáln´ı optimalizace. Je to moˇzno ˇreˇsit pouˇzit´ım Grayova kódu [12], kdy ˇretˇezce reprezentuj´ıc´ı sousedn´ı celoˇc´ıselné hodnoty se liˇs´ı jen v jednom bitu, ale pˇrináˇs´ı to dalˇs´ı ˇcasov´ y nárok na konverzi do Grayova kódu a zpˇet. Byly navrˇzeny i varianty genetick´ ych algoritm˚ u s reprezentac´ı jedince jako vektoru ˇc´ıseln´ ych hodnot v pohyblivé ˇcárce. Pak je ovˇsem nutné uˇz´ıvat jiné operace kˇr´ıˇzen´ı a mutace neˇz ty, které se pouˇz´ıvaj´ı pro bitové ˇretˇezce. Nˇekteré takové operace si ukáˇzeme na dalˇs´ıch evoluˇcn´ıch algoritmech v tomto textu. Shrnut´ı: - Evoluce v biologick´ ych systémech - Genetická informace, chromozóm, fitness - Operátory selekce, kˇr´ıˇzen´ı a mutace v genetick´ ych algoritmech

Kontroln´ı ot´ azky: 1. Proˇc pˇri hledán´ı globáln´ıho minima nelze uˇz´ıvat jako fitness hodnotu u ´ˇcelové funkce? ˇ ım se liˇs´ı dvoubodové kˇr´ıˇzen´ı od jednobodového? 2. C´ 3. Jaká je role mutace v genetick´ ych algoritmech?


36

6.2

Evoluˇ cn´ı strategie

Pr˚ uvodce studiem: V této kapitole jsou vysvˇetleny základy evoluˇcn´ı strategie, která patˇr´ı k nejstarˇs´ım a stále velmi populárn´ım a inspiruj´ıc´ım evoluˇcn´ım algoritm˚ um. Poˇc´ıtejte asi se tˇremi aˇz ˇctyˇrmi hodinami studia.

P˚ uvodn´ı nejjednoduˇsˇs´ı verzi algoritmu navrhli v ˇsedesát´ ych létech minulého stolet´ı Schwefel a Rechenberg, kdyˇz potˇrebovali nalézt optimáln´ı tvar obtékaného tˇelesa. Základn´ı myˇslenka je podobná slepému náhodnému prohledáván´ı, ale rozd´ıl je v tom, ˇze nov´ y bod y se generuje jako mutace bodu x tak, ˇze jednotlivé sloˇzky vektoru x se zmˇen´ı pˇriˇcten´ım hodnot normálnˇe rozdˇelen´ ych náhodn´ ych veliˇcin y = x + u,

u ∼ N (0, σ 2 I),

(6.1)

tj. u = (U1 , U2 , . . . , Ud ) je náhodn´ y vektor, jehoˇz kaˇzd´ y prvek je náhodná veliˇcina 2 Ui ∼ N (0, σ ), i = 1, 2, . . . , d a tyto veliˇciny jsou navzájem nezávislé. Algoritmus lze zapsat v Matlabu takto: function [x,fx] = es1p1(fn_name,a,b,t,sigma) % input parameters: % fn_name function to be minimized (M file) % a, b row vectors, limits of search space % t number of iterations % sigma standard deviation (const) % output: % fx the minimal function value found % x minimum found by search d = length(a); fx = realmax; x = a + (b - a).*rand(1,d); for i = 1:t y = x + sigma*randn(1,d); y = zrcad(y, a, b); fy = feval(fn_name, y); if fy < fx x = y; fx = fy; end end

6.2 Evoluˇcn´ı strategie

37

Porovnejte tento zdrojov´ y text se zdrojov´ ym textem algoritmu slepého prohledáván´ı v kap. 4 a uvˇedomte si rozd´ıly. Evoluˇcn´ı terminologi´ı m˚ uˇze b´ yt tento algoritmus popsán tak, ˇze z rodiˇcovské generace velikosti 1 vzniká generace potomk˚ u rovnˇeˇz velikosti 1, potomek vzniká mutac´ı z jednoho rodiˇce podle rov. (6.1) a operátorem selekce je v´ ybˇer lepˇs´ıho jedince z dvojice (tzv. turnajov´ y v´ ybˇer). Tento algoritmus je v literatuˇre oznaˇcován zkratkou ES(1+1). Pˇ r´ıklad 6.1 Porovnáme u ´ˇcinnost této nejjednoduˇsˇs´ı varianty evoluˇcn´ı strategie se slep´ ym náhodn´ ym prohledáván´ım. Jak ze srovnán´ı zdrojového kódu obou algoritm˚ u vid´ıme, jsou velmi podobné, evoluˇcn´ı strategie má o jeden vstupn´ı parametr v´ıc, a to o smˇerodatnou odchylku sigma, která je shodná pro vˇsechny sloˇzky nového generovaného bodu y a konstantn´ı pro celé prohledáván´ı. V implementaci funkce es1p1 je nav´ıc uˇzito zrcadlen´ı, aby novˇe vygenerovan´ y bod nebyl vnˇe prostoru D. Pro spouˇstˇen´ı experimentáln´ıho porovnán´ı si nap´ıˇseme následuj´ıc´ı skript, kter´ y uspoˇr´ı pracné opakované zadáván´ı pˇr´ıkaz˚ u z klávesnice a umoˇzn´ı i zaznamenat dosaˇzené v´ ysledky z jednotliv´ ych bˇeh˚ u do souboru, kter´ y pak bude vstupem pro statistické zpracován´ı v´ ysledk˚ u. V´ ybˇer u ´lohy k ˇreˇsen´ı zaˇr´ıd´ıme vˇzdy odstranˇen´ım znaku % na zaˇca´tku zvoleného ˇra´dku a “vyprocentován´ım” ostatn´ıch ˇra´dk˚ u se jménem funkce a dalˇs´ımi vstupn´ımi u ´daji. Podobn´ y skript si m˚ uˇzete snadno vytvoˇrit i pro spouˇstˇen´ı slepého prohledáván´ı. Testy provedeme na ˇsesti testovac´ıch funkc´ıch uveden´ ych v kap. 3.2, dimenze problému d = 2, funkce jsou uˇzity bez posunu (o = (0, 0), viz kap. 3.2), pro vˇsechny funkce je prohledáván´ı ukonˇceno po 10 000 vyhodnocen´ı funkce a pro kaˇzdou u ´lohu se provede 10 opakován´ı. Volba hodnoty sigma je zˇrejmá ze zdrojového kódu spouˇstˇec´ıho skriptu. % start tol=1e-4; % tolerance f(x) pro f_near cfn=0; s=420.968746; % pro schwefel0 % fn_name=’ackley’; b=30*ones(1,2); cxst=zeros(1,2); f_near=cfn+tol; % fn_name=’ackley’; b=30*ones(1,5); cxst=zeros(1,5); f_near=cfn+tol; % fn_name=’ackley’; b=30*ones(1,10); cxst=zeros(1,10); f_near=cfn+tol; % fn_name=’ackley’; b=30*ones(1,30); cxst=zeros(1,30); f_near=cfn+tol; fn_name=’rosen’; b=[2.048 2.048]; cxst=ones(1,2); f_near=cfn+tol; % fn_name=’rosen’; b=2.048*ones(1,5); cxst=ones(1,5);f_near=cfn+tol; % fn_name=’rosen’; b=2.048*ones(1,10); cxst=ones(1,10);f_near=cfn+tol; % fn_name=’rosen’; b=2.048*ones(1,30); cxst=ones(1,30);f_near=cfn+tol; % fn_name=’dejong1’; b=5.12*ones(1,2); cxst=zeros(1,2);f_near=cfn+tol;

38


% fn_name=’dejong1’; b=5.12*ones(1,5); cxst=zeros(1,5);f_near=cfn+tol; % fn_name=’dejong1’; b=5.12*ones(1,10);cxst=zeros(1,10);f_near=cfn+tol; % fn_name=’dejong1’; b=5.12*ones(1,30);cxst=zeros(1,30);f_near=cfn+tol; % fn_name=’griewank’; b=400*ones(1,2); cxst=zeros(1,2);f_near=cfn+tol; % fn_name=’griewank’; b=400*ones(1,5); cxst=zeros(1,5);f_near=cfn+tol; % fn_name=’griewank’; b=400*ones(1,10); cxst=zeros(1,10);f_near=cfn+tol; % fn_name=’griewank’; b=400*ones(1,30); cxst=zeros(1,30);f_near=cfn+tol; % fn_name=’schwefel0’;b=500*ones(1,2); cxst=s*ones(1,2);f_near=cfn+tol; % fn_name=’schwefel0’;b=500*ones(1,5); cxst=s*ones(1,5);f_near=cfn+tol; % fn_name=’schwefel0’;b=500*ones(1,10);cxst=s*ones(1,10);f_near=cfn+tol; % fn_name=’schwefel0’;b=500*ones(1,30);cxst=s*ones(1,30);f_near=cfn+tol; % fn_name=’rastrig’; b=5.12*ones(1,2); cxst=zeros(1,2);f_near=cfn+tol; % fn_name=’rastrig’; b=5.12*ones(1,5); cxst=zeros(1,5);f_near=cfn+tol; % fn_name=’rastrig’; b=5.12*ones(1,10);cxst=zeros(1,10);f_near=cfn+tol; % fn_name=’rastrig’; b=5.12*ones(1,30);cxst=zeros(1,30);f_near=cfn+tol; max_evals=5000; a=-b; d=length(a); evals= d * max_evals; sigma=(b - a) / 10; sigma = mean(sigma); fid=fopen(’es11.002’,’a’); disp(fn_name); for krok=1:10 disp(’krok’); disp(krok); [xmin,fnmin] = es1p1(fn_name, a, b, evals, sigma); fprintf(fid,’%-15s’, ’ES11’); fprintf(fid,’%-10s %4.0f’, fn_name, d); fprintf(fid, ’%4.0f ’, krok); fprintf(fid, ’%8.4e ’, sigma); fprintf(fid,’%8.0f’, evals); fprintf(fid,’ %22.11e’, fnmin); fprintf(fid,’%1s\n’,’ ’); end fclose(fid);


39

V kaˇzdé u ´loze vyhodnocujeme nalezenou pr˚ umˇernou hodnotu funkce spoˇc´ıtanou z 10 opakován´ı. V´ ysledky jsou uvedeny v následuj´ıc´ı tabulce.

blindsearch fprum

ES(1+1)

std

fprum

std

ackley

0.6214 0.2528

0.2719

0.1937

dejong1

0.0069 0.0080

0.0002

0.0001

griewank

0.1155 0.0497

0.0291

0.0213

rastrig

0.3306 0.3900

0.0531

0.0654

rosen

0.0049 0.0083

0.0004

0.0004

schwefel0 4.3406 3.3183 130.5089 87.5699

Vid´ıme, ˇze s v´ yjimkou Schwefelovy funkce dosáhla i tato jednoduchá varianta evoluˇcn´ı strategie v´ yraznˇe lepˇs´ıch v´ ysledk˚ u. Pr˚ umˇer nalezen´ ych funkˇcn´ıch hodnot je v´ yraznˇe menˇs´ı a také smˇerodatná odchylka nalezen´ ych funkˇcn´ıch hodnot je podstatnˇe menˇs´ı neˇz u slepého prohledáván´ı. Horˇs´ı v´ ysledky u Schwefelovy funkce jsou zp˚ usobeny t´ım, ˇze prohledáván´ı konˇc´ı v oblasti nˇekterého z lokáln´ıch minim a tuto oblast nen´ı schopno opustit. Algoritmus ES(1+1), kter´ y vyuˇz´ıvá své “znalosti” z´ıskané bˇehem hledán´ı (byt’ jen velmi omezenˇe), je u ´ˇcinnˇejˇs´ı neˇz u ´plnˇe slepé a bezhlavé prohledáván´ı. Konec pˇ r´ıkladu.

Zaj´ımavou otázkou je to, jak volit hodnoty smˇerodatn´ ych odchylek pro mutaci. V rov. (6.1) jsou vˇsechny hodnoty σj , j = 1, 2, . . . , d shodné a konstantn´ı po celé vyhledáván´ı. To samozˇrejmˇe nen´ı nutné, kaˇzdá dimenze m˚ uˇze m´ıt svou hodnotu smˇerodatné odchylky, tedy budeme pracovat s vektorem smˇerodatn´ ych odchylek σ = (σ1 , σ2 , . . . , σd ) Hodnoty smˇerodatn´ ych odchylek mohou b´ yt urˇceny jako pomˇerná ˇcást velikosti hrany d-rozmˇerného kvádru D, tj. σ = c(b − a), kde 0 < c < 1 je zadaná konstanta. Je zˇrejmé, ˇze ˇc´ım je hodnota c vˇetˇs´ı, t´ım d˚ ukladnˇeji bude doména D prohledávána, ale za cenu pomalejˇs´ı konvergence, naopak pˇr´ıliˇs malá hodnota c zv´ yˇs´ı riziko pˇredˇcasné konvergence algoritmu v lokáln´ım minimu.


40

Nav´ıc, vektor smˇerodatn´ ych odchylek nemus´ı b´ yt v kaˇzdém iteraˇcn´ım kroku stejn´ y, ale hodnoty jeho prvk˚ u se mohou adaptovat podle pr˚ ubˇehu procesu vyhledáván´ı. Rechenberg studoval vliv velikosti mutace pˇri generován´ı potomka, tj. hodnot σ na rychlost konvergence algoritmu a na dvou jednoduch´ ych funkc´ıch odvodil tzv. pravidlo jedné pˇetiny u ´spˇeˇsnosti. Empiricky pak byla ovˇeˇrena uˇziteˇcnost tohoto pravidla i pro jiné funkce. Hodnoty smˇerodatn´ ych odchylek se v i-té iteraci hledán´ı upravuj´ı podle pravidla vyjádˇreného rovnic´ı  1    c1 σi−1 kdyˇz φ(n) < 5 σi = (6.2) c2 σi−1 kdyˇz φ(n) > 51    σ jinak i−1

c1 < 1 a c2 > 1 jsou vstupn´ı parametry, kter´ ymi se ˇr´ıd´ı zmenˇsován´ı ˇci zvˇetˇsován´ı hodnot smˇerodatn´ ych odchylek podle relativn´ı ˇcetnosti u ´spˇechu φ(n) v pˇredcházej´ı´ echem se rozum´ı to, ˇze f (y) < f (x). Poˇcet krok˚ c´ıch n kroc´ıch. Uspˇ u n je vstupn´ı . parametr. Obvykle se vol´ı hodnoty c1 = 0.82 a c2 = 1/c1 = 1.22. Pod vlivem rozvoje jin´ ych evoluˇcn´ıch algoritm˚ u bylo dalˇs´ım krokem v evoluˇcn´ı strategii zaveden´ı populace velikosti vˇetˇs´ı neˇz jedna. Oznaˇc´ıme-li velikost rodiˇcovské populace µ a velikost populace potomk˚ u λ, λ > µ, pak m˚ uˇzeme uvaˇzovat o dvou variantách algoritmu evoluˇcn´ı strategie, v literatuˇre obvykle oznaˇcovan´ ych jako ES(µ+λ) a ES(µ, λ): • ES(µ + λ) - generace potomk˚ u je vytvoˇrena µ jedinci s nejmenˇs´ımi funkˇcn´ımi hodnotami ze vˇsech µ + λ jedinc˚ u jak rodiˇcovské populace, tak populace potomk˚ u. V této variantˇe se dodrˇzuje tzv. elitismus, tj. nejlepˇs´ı dosud nalezen´ y jedinec vˇzdy pˇreˇz´ıvá a pˇrecház´ı do nové generace. • ES(µ, λ) - generace potomk˚ u je vytvoˇrena µ jedinci s nejmenˇs´ımi funkˇcn´ımi hodnotami z λ jedinc˚ u populace potomk˚ u. Rodiˇcovská generace je tedy kompletnˇe nahrazena nejlepˇs´ımi jedinci z populace potomk˚ u. Doporuˇcuje se, aby velikost populace potomk˚ u λ byla volena nˇekolikanásobnˇe vˇetˇs´ı neˇz velikost rodiˇcovské populace µ. Také se uvád´ı, ˇze varianta ES(µ, λ) obvykle konverguje pomaleji neˇz ES(µ+λ), ale s menˇs´ı tendenc´ı ukonˇcit prohledáván´ı v lokáln´ım minimu. Implementace algoritmu ES(µ+λ) s adaptac´ı smˇerodatn´ ych odchylek podle pravidla jedné pˇetiny v Matlabu je uvedena v následuj´ıc´ım v´ ypisu: function [x, fx, func_evals] =... es_mpl(fn_name, a, b,... M, k, lfront, my_eps, max_evals, shift)


d = length(a); sigma = (b - a) / 6; % initial values of sigma fr = rand(1, lfront) < 0.2; % initial values in last lfront steps % fr je vektor delky lfront % s prvky 1 nebo 0 indikujicimi uspech nebo neuspech v poslednich % lfront generovani bodu y L = M*k; % k>1, integer, velikost generace potomku P = zeros(M, d+1); % generace rodicu Q = zeros(L, d+1); % generace potomku for i=1:M P(i,1:d) = a + (b-a).*rand(1,d); P(i,d+1)= feval(fn_name, P(i,1:d)- shift); end % 0-th generation initialized P = sortrows(P, d+1); func_evals = M; while ((P(M,d+1) - P(1,d+1)) > my_eps)&&... (func_evals <= d*max_evals - L) % main loop for t = 1:k for i = 1:M x = P(i, 1:d); y = x + sigma.*randn(1,d); y = zrcad(y,a,b); fy = feval(fn_name, y - shift); func_evals = func_evals + 1; Q((t-1)*M + i, :) = [y fy]; fr(1)= []; if fy < P(i,d+1) % indikator uspechu a neuspechu fr(end + 1) = 1; else fr(end + 1) = 0; end end % adaptace sigma podle cetnosti uspechu a neuspechu % v poslednich lfront krocich if sum(fr) / lfront < 0.2 sigma = sigma * 0.82; elseif sum(fr) / lfront > 0.2 sigma = sigma * 1.22; end end % vytvorena generace potomku velikosti L

41


42

PQ =[ P; Q]; % spojeni generace rodicu a potomku PQ = sortrows(PQ, d+1); P = PQ(1:M,:); % nova rodicovska generace end % main loop - end x = P(1,1:d); fx = P(1,d+1); Autoˇri algoritmu evoluˇcn´ı strategie v minulosti museli ˇcelit v´ ytkám, ˇze evoluˇcn´ı strategie z osvˇedˇcen´ ych evoluˇcn´ıch operátor˚ u v˚ ubec nevyuˇz´ıvá kˇr´ıˇzen´ı. Proto byly navrˇzeny pokroˇcilejˇs´ı varianty evoluˇcn´ı strategie, ve kter´ ych je kˇr´ıˇzen´ı obsaˇzeno. Základn´ı idea takového kˇr´ıˇzen´ı spoˇc´ıvá v tom, ˇze u kaˇzdého jedince v populaci je kromˇe souˇradnic v prostoru D uchováván i jeho vektor smˇerodatn´ ych odchylek. Vektor smˇerodatn´ ych odchylek potomka se pak generuje kˇr´ıˇzen´ım s náhodnˇe vybran´ ym jin´ ym r s jedincem. Máme-li dva rodiˇce r, s s vektory smˇerodatn´ ych odchylek σ , σ , pak vektor smˇerodatn´ ych odchylek jejich potomka je urˇcován napˇr. jako σ=

σr + σs , 2

coˇz je tzv. kˇr´ıˇzen´ı pr˚ umˇerem, nebo podle nˇejakého podobného pravidla pro kˇr´ıˇzen´ı. Pro funkce, u kter´ ych se dá oˇcekávat globáln´ı minimum v protáhlém u ´dol´ı, jehoˇz smˇer nen´ı rovnobˇeˇzn´ y se ˇza´dnou dimenz´ı prostoru D, je vhodné uˇz´ıt jeˇstˇe sofistikovanˇejˇs´ı variantu evoluˇcn´ı strategie, v n´ıˇz u kaˇzdého jedince v populaci se kromˇe souˇradnic v prostoru D a jeho vektoru smˇerodatn´ ych odchylek uchovávaj´ı i mimodiagonáln´ı prvky kovarianˇcn´ı matice rozmˇeru (d × d). Nov´ y bod vzniká mutac´ı, tj. pˇriˇcten´ım náhodného vektoru u ∼ N (0, Σ), kde Σ je kovarianˇcn´ı matice vektor˚ u souˇradnic bod˚ u v populaci. Také kˇr´ıˇzen´ı m˚ uˇze prob´ıhat nejen na vektorech smˇerodatn´ ych odchylek, ale na celé kovarianˇcn´ı matici. Takové modifikace evoluˇcn´ı strategie jsou v literatuˇre oznaˇcovány zkratkou CMES (Covariance Matrix Evolution Strategy). Jsou vˇsak nejenom implementaˇcnˇe podstatnˇe nároˇcnˇejˇs´ı, ale také maj´ı vˇetˇs´ı pamˇet’ové nároky i vˇetˇs´ı ˇcasovou spotˇrebu na jeden iteraˇcn´ı krok. Proto se uˇz´ıvaj´ı sp´ıˇse jen u problém˚ u, kdy vyhodnocen´ı u ´ˇcelové funkce je ˇcasovˇe nároˇcné a sn´ıˇzen´ı poˇctu vyhodnocen´ı funkce uspoˇr´ı v´ıce ˇcasu, neˇz kter´ y je spotˇrebován na opakované vyhodnocován´ı kovarianˇcn´ı matice a sloˇzitˇejˇs´ı kˇr´ıˇzen´ı. M´ısto kovarianc´ı dvojic souˇradnic se mohou uˇz´ıvat smˇerové u ´hly, které lze z kovarianˇcn´ı matice snadno vyhodnotit. Porovnáme-li evoluˇcn´ı strategii ES (µ + λ) s algoritmem CRS, mohli bychom CRS oznaˇcit jako evoluˇcn´ı algoritmus typu (µ + 1), nebot’ nová generace má jen jednoho jedince, kter´ y m˚ uˇze nahradit nejhorˇs´ıho jedince v populaci. Stejnˇe jako v ES(µ + λ) se spojuje stará a nová generace a z této spojené populace pˇreˇz´ıvá µ nejlepˇs´ıch jedinc˚ u, tedy uˇz´ıvá se elitismus. Na rozd´ıl od evoluˇcn´ı strategie se nov´ y bod v CRS negeneruje mutac´ı pˇriˇcten´ım normálnˇe rozdˇeleného náhodného vektoru, ale lokáln´ı


43

heuristikou, napˇr. reflex´ı náhodnˇe vybraného simplexu, kterou m˚ uˇzeme povaˇzovat za kombinaci mutace a kˇr´ıˇzen´ı. Shrnut´ı: - mutace v evoluˇcn´ı strategii, jej´ı ˇr´ıdic´ı parametry - algoritmus ES (1+1) - adaptace hodnot smˇerodatn´ ych odchylek - ES(µ + λ) a ES(µ, λ), elitismus

Kontroln´ı ot´ azky: 1. Co je pravidlo jedné pˇetiny a k ˇcemu se uˇz´ıvá? 2. Porovnejte v´ yhody a nev´ yhody variant ES(µ, λ) a ES(µ + λ). 3. Jak je pravidlo jedné pˇetiny implementováno ve variantˇe ES(µ + λ) v Matlabu a kter´ ym parametrem ˇr´ıd´ıme délku sledovaného poˇctu pˇredchoz´ıch iterac´ı? 4. Jak byste porovnali u ´ˇcinnosti algoritmu ES(µ+λ) a ES (1+1) experimentálnˇe? O kterém z tˇechto algoritm˚ u se domn´ıváte, ˇze bude u ´ˇcinnˇejˇs´ı?


44

6.3

Diferenci´ aln´ı evoluce

Pr˚ uvodce studiem: V této kapitole se seznám´ıte s algoritmem diferenciáln´ı evoluce. Tento algoritmus byl navrˇzen nedávno a poprvé publikován v roce 1995. Bˇehem nˇekolika let se stal velmi populárn´ı, byl podrobnˇe studován na mnoha pracoviˇst´ıch a je ˇcasto aplikován na problémy hledán´ı globáln´ıho minima. Je to pˇr´ıklad jednoduchého heuristického hledán´ı, ve kterém se uˇz´ıvaj´ı evoluˇcn´ı operátory. Poˇc´ıtejte asi se tˇremi aˇz pˇeti hodinami studia.

Diferenciáln´ı evoluce (DE) je postup heuristického hledán´ı minima funkc´ı, kter´ y ych letech. Experimentáln´ı v´ ysledky navrhli R. Storn a K. Price [20] v devadesát´ z testován´ı i zkuˇsenosti z ˇcetn´ ych aplikac´ı ukazuj´ı, ˇze ˇcasto konverguje rychleji neˇz jiné stochastické algoritmy pro globáln´ı optimalizaci. Algoritmus diferenciáln´ı evoluce vytváˇr´ı novou generaci Q tak, ˇze postupnˇe pro kaˇzd´ y bod xi ze staré generace P vytvoˇr´ı jeho potenciáln´ıho konkurenta y a do nové populace z této dvojice zaˇrad´ı bod s niˇzˇs´ı funkˇcn´ı hodnotou. Algoritmus m˚ uˇzeme v pseudokódu zapsat takto:

generuj poˇcáteˇcn´ı populaci P (N bod˚ u náhodnˇe v D) repeat for i := 1 to N do generuj vektor u mutac´ı vytvoˇr vektor y kˇr´ıˇzen´ım u a xi if f (y) < f (xi ) then do Q zaˇrad’ y else do Q zaˇrad’ xi endfor P := Q until podm´ınka ukonˇcen´ı

Generován´ı vektoru u pˇredstavuje v DE mutaci, tento vektor-mutant je pak jedn´ım z rodiˇcovsk´ ych vektor˚ u pro kˇr´ıˇzen´ı, druh´ ym rodiˇcem je aktuáln´ı protˇejˇsek ve staré generaci, nov´ y zkusm´ y bod y je tedy produktem evoluˇcn´ıch operátor˚ u mutace a kˇr´ıˇzen´ı, selekce se provád´ı turnajem s jeho protˇejˇskem ve staré generaci. Je nˇekolik moˇzn´ ych zp˚ usob˚ u mutace a dva druhy kˇr´ıˇzen´ı. Varianty tˇechto r˚ uzn´ ych strategi´ı diferenciáln´ı evoluce se ˇcasto oznaˇcuj´ı zkratkou DE/m/a/k, kde m znamená uˇzit´ y zp˚ usob mutace, a je poˇcet rozd´ıl˚ u (diferenc´ı) vektor˚ u v operaci mutace a k je typ kˇr´ıˇzen´ı.

6.3 Diferenciáln´ı evoluce

45

Nejˇcastˇeji uˇz´ıvanou mutac´ı v DE je postup oznaˇcovan´ y rand/1/, kter´ y generuje bod u ze tˇr´ı bod˚ u ze staré populace podle vztahu u = r 1 + F (r 2 − r 3 ),

(6.1)

r 1 , r 2 , r 3 jsou navzájem r˚ uzné body náhodnˇe vybrané z populace P r˚ uzné od aktuáln´ıho bodu xi , F > 0 je vstupn´ı parametr. Generován´ı bodu u podle rovnice (6.1) je graficky znázornˇeno následuj´ıc´ım obrázkem.

10 9

D

8

F(r −r ) 2 3

7

r

u

1

6 5

xi

r2−r3

r3

4

r

2

3 2 1 0 −1

0

2

4

6

8

10

Dalˇs´ı ˇcasto uˇz´ıvan´ y druh mutace vyuˇz´ıvá nejlepˇs´ı bod z populace P a ˇctyˇri dalˇs´ı náhodnˇe vybrané body podle vztahu u = xbest + F (r 1 + r 2 − r 3 − r 4 ),

(6.2)

kde r 1 , r 2 , r 3 , r 4 jsou navzájem r˚ uzné body náhodnˇe vybrané z populace P r˚ uzné od aktuáln´ıho bodu xi i od bodu xbest s nejniˇzˇs´ı funkˇcn´ı hodnotou v populaci P . F > 0 je opˇet vstupn´ı parametr. V (6.2) vid´ıme, ˇze k vektoru xbest se pˇriˇc´ıtaj´ı dva rozd´ıly náhodnˇe vybran´ ych vektor˚ u násobené F , proto je tato mutace oznaˇcována best/2/. Kaelo a Ali [11] navrhli modifikovat mutaci rand/1/ tak, ˇze vektor r 1 v (6.1) je ten s nejmenˇs´ı funkˇcn´ı hodnotou mezi vektory r 1 , r 2 a r 3 . Oznaˇcuj´ı takovou mutaci jako náhodnou lokalizaci (random localization) a tuto mutaci m˚ uˇzeme oznaˇcit jako randrl/1/. V rozsáhl´ ych experimentáln´ıch testech ukázali, ˇze oproti rand/1/ tato mutace sniˇzuje poˇcet vyhodnocen´ı potˇrebn´ ych k dosaˇzen´ı podm´ınky ukonˇcen´ı aˇz o tˇretinu a vˇetˇsinou je zachována spolehlivost dobrého pˇribl´ıˇzen´ı ke globáln´ımu minimu.

46


V posledn´ı dobˇe se v nˇekter´ ych variantách DE uˇz´ıvá mutace current-to-best/2/, která generuje vektor u podle vztahu u = xi + F (xbest − xi ) + F (r 1 − r 2 )

(6.3)

Tato mutace má podobnˇe jako (6.2) tendenci zrychlovat konvergenci algoritmu a zvyˇsovat riziko pˇredˇcasné konvergence v lokáln´ım minimu. Kromˇe uveden´ ych druh˚ u mutace se v literatuˇre objevuje i celá ˇrada dalˇs´ıch, ale mutace rand/1/ je stále v aplikac´ıch DE nejˇcastˇeji uˇz´ıvanou. Nov´ y vektor y vznikne kˇr´ıˇzen´ım vektoru u a vektoru xi . V diferenciáln´ı evoluci se uˇz´ıvaj´ı dva typy kˇr´ıˇzen´ı, a to binomické kˇr´ıˇzen´ı a exponenci´ aln´ı kˇr´ıˇzen´ı. Oba typy kˇr´ıˇzen´ı byly navrˇzeny v prvn´ıch publikac´ıch o DE [20]. V binomickém kˇr´ıˇzen´ı vznikne vektor y tak, ˇze kter´ ykoli prvek tohoto vektoru m˚ uˇze roven bud’ odpov´ıdaj´ıc´ımu prvku vektoru xi nebo odpov´ıdaj´ıc´ımu prvku vektoru u. Pravdˇepodobnost, ˇze ve vektoru y bude prvek vektoru u je u ´mˇerná hodnotˇe zadaného vstupn´ıho parametru CR. Binomické kˇr´ıˇzen´ı prob´ıhá podle následuj´ıc´ıho pravidla: { uj kdyˇz Rj ≤ CR nebo j = I yj = (6.4) xij kdyˇz Rj > CR a j ̸= I , kde I je náhodnˇe vybrané celé celé ˇc´ıslo z {1, 2, . . . , d}, Rj ∈ (0, 1) jsou voleny náhodnˇe a nezávisle pro kaˇzdé j a CR ∈ [0, 1] je vstupn´ı parametr. Z pravidla (6.4) vid´ıme, ˇze nejménˇe jeden prvek vektoru u pˇrecház´ı do nového bodu zkusmého y, dokonce i pˇri volbˇe CR = 0. Strategie diferenciáln´ı evoluce uˇz´ıvaj´ıc´ı binomiáln´ı kˇr´ıˇzen´ı se oznaˇcuj´ı zkratkou DE/m/a/bin. V exponenciáln´ım kˇr´ıˇzen´ı se vymˇen ˇuje náhodnˇe urˇcen´ y poˇcet sousedn´ıch prvk˚ u vektor˚ u. V tom je exponenciáln´ı kˇr´ıˇzen´ı podobné dvoubodovému kˇr´ıˇzen´ı popsaném v genetick´ ych algoritmech. Poˇcáteˇcn´ı pozice pro kˇr´ıˇzen´ı se vybere náhodnˇe z {1, . . . , d} a pak L následuj´ıc´ıch prvk˚ u (poˇc´ıtáno cyklicky, tj. za d-t´ ym prvkem následuje prvn´ı prvek vektoru) se vloˇz´ı z vektoru u. Jeden prvek vektoru u se vloˇz´ı vˇzdy, pravdˇepodobnost vloˇzen´ı dalˇs´ıch prvk˚ u pak klesá exponenciálnˇe. Strategie diferenciáln´ı evoluce uˇz´ıvaj´ıc´ı exponenciáln´ı kˇr´ıˇzen´ı se oznaˇcuj´ı zkratkou DE/m/a/exp. Popsané exponenciáln´ı kˇr´ıˇzen´ı lze snadno implementovat vyuˇzit´ım cyklu while pro zvˇetˇsován´ı poˇctu prvk˚ u z vektoru u, jak je ukázáno v následuj´ıc´ım zdrojovém textu funkce rand1_exp, která vytvoˇr´ı nov´ y pokusn´ y bod y strategi´ı DE/rand/1/exp.


47

% rand/1/exp function y = rand1_exp(P, F, CR, expt) N = length(P(:, 1)); d = length(P(1,:)) - 1; y = P(expt(1),1:d); % expt(1) je index akt. radku v P vyb = nahvyb_expt(N, 3, expt); % three random points without expt r1 = P(vyb(1),1:d); r2 = P(vyb(2),1:d); r3 = P(vyb(3),1:d); u = r1 + F*(r2 - r3); nah_pozice = 1 + fix(d*rand(1)); pridej = 0; while rand(1) < CR && pridej < d - 1 pridej = pridej + 1; end delka_nazacatku = nah_pozice + pridej - d; if delka_nazacatku <= 0 indexy_zmeny = nah_pozice: nah_pozice + pridej; else indexy_zmeny = [1: delka_nazacatku, nah_pozice:d]; end y(indexy_zmeny)= u(indexy_zmeny);

Pravdˇepodobnost pm zaˇrazen´ı prvk˚ u mutaˇcn´ıho vektoru u do zkusmého vektoru y je moˇzno definovat jako stˇredn´ı hodnotu relativn´ı ˇcetnosti tˇechto prvk˚ u, tj. pm = E(L)/d. Vztah mezi touto pravdˇepodobnost´ı pm a ˇr´ıdic´ım parametrem CR podrobnˇe studovala Zaharie [28]. U binomického kˇr´ıˇzen´ı je tato závislost lineárn´ı, pm = CR(1 − 1/d) + 1/d,

(6.5)

zat´ımco pro exponenciáln´ı kˇr´ıˇzen´ı je tato závislost silnˇe nelineárn´ı a odchylka od linearity se zvˇetˇsuje s rostouc´ı dimenz´ı prohledáváného prostoru. pm =

1 − CRd , d(1 − CR)

pro CR < 1.

(6.6)

Graf závislosti CR a pm pro d = 18 je na následuj´ıc´ım obrázku. Nelinearitu této závislosti pro exponenciáln´ı kˇr´ıˇzen´ı je nutno vz´ıt v u ´vahu pˇri zadáván´ı hodnoty . parametru CR. Napˇr. pro tuto dimenzi problému, chceme-li dosáhnout pm = 0.5, je . nutno zadat CR = 0.9.


48

1 Exponential 0.8

0.6 CR

Binomial 0.4

0.2

p

1

0 0

0.2

p2

0.4 0.6 Probability of mutation

p3

0.8

1

V´ yhodou algoritmu diferenciáln´ı evoluce je jeho jednoduchost a v´ ypoˇcetnˇe nenároˇcné generován´ı nového bodu y, které lze nav´ıc velmi efektivnˇe implementovat. Pro zvolenou strategii DE/m/a/k je potˇreba zadat jen tˇri hodnoty ˇr´ıdic´ıch parametr˚ u, a to velikost populace N , hodnotu parametru F pro mutaci a hodnotu parametru CR pro kˇr´ıˇzen´ı. Nev´ yhodou je pomˇernˇe velká citlivost algoritmu na nastaven´ı hodnot parametr˚ u F a CR. Storn a Price doporuˇcuj´ı volit hodnoty N = 10d, F = 0.8 a CR = 0.5 a pak tyto hodnoty modifikovat podle empirické zkuˇsenosti z pozorovaného pr˚ ubˇehu hledán´ı. To je ovˇsem doporuˇcen´ı velmi vágn´ı a ˇcasovˇe nároˇcné. D˚ uleˇzitá je i zkuˇsenost a dobrá intuice ˇreˇsitele problému. Vˇetˇsinou lze uˇz´ıt hodnoty 0.5 ≤ F ≤ 1 a hodnoty CR z celého oboru, 0 ≤ CR ≤ 1. Také velikost populace N ˇcasto postaˇc´ı menˇs´ı neˇz N = 10d, t´ım je moˇzné zv´ yˇsit rychlost konvergence, ale souˇcasnˇe se pˇri menˇs´ı velikosti populace zvyˇsuje riziko pˇredˇcasné konvergence. Citlivost algoritmu DE na nastaven´ı hodnot ˇr´ıdic´ıch parametr˚ u vedla k r˚ uzn´ ym návrh˚ um adaptivn´ıch verz´ı DE, ve kter´ ych se hodnoty ˇr´ıdic´ıch parametr˚ u mˇen´ı adaptivnˇe podle ˇreˇseného problému v pr˚ ubˇehu prohledáván´ı, takˇze uˇzivatel pouze nastav´ı jejich poˇca´teˇcn´ı hodnoty nebo dokonce se nemus´ı o nastaven´ı parametr˚ u starat v˚ ubec. Pˇr´ıkladem takového postupu je adaptivn´ı nastaven´ı parametru F , kter´ y navrhli Ali a Törn. Hodnota parametru F se v kaˇzdé generaci mˇen´ı podle adaptivn´ıho pravidla { max(Fmin , 1 − | ffmax |) if | ffmax |<1 min min F = (6.7) min max(Fmin , 1 − | ffmax |) jinak, kde fmin , fmax jsou minimáln´ı a maximáln´ı funkˇcn´ı hodnoty v populaci P a Fmin je vstupn´ı parametr, kter´ y zabezpeˇcuje, aby bylo F ∈ [Fmin , 1). Autoˇri doporuˇcuj´ı volit hodnotu Fmin ≥ 0.45. Pˇredpokládá se, ˇze tento zp˚ usob v´ ypoˇctu F udrˇzuje prohledáván´ı diverzifikované v poˇca´teˇcn´ım stádiu a intenzivnˇejˇs´ı v pozdˇejˇs´ı fázi prohledáván´ı, coˇz obvykle zvyˇsuje spolehlivost hledán´ı i rychlost konvergence. Nˇekteré


49

dalˇs´ı moˇznosti adaptace vyhledávac´ı strategie pro hledán´ı globáln´ıho minima ukáˇzeme v závˇereˇcné kapitole tohoto uˇcebn´ıho textu. Porovnáme-li diferenciáln´ı evoluci s algoritmem CRS, vid´ıme, ˇze v diferenciáln´ı evoluci se nenahrazuje nejhorˇs´ı bod v populaci, ale pouze horˇs´ı bod ve dvojici, coˇz zabezpeˇcuje vˇetˇs´ı diverzitu bod˚ u populace po delˇs´ı obdob´ı, takˇze pˇri stejné velikosti populace DE ve srovnán´ı s CRS má vˇetˇsinou menˇs´ı tendenci ukonˇcit prohledáván´ı v lokáln´ım minimu, ale za to plat´ıme pomalejˇs´ı konvergenc´ı pˇri stejné podm´ınce ukonˇcen´ı. Nejˇcastˇeji uˇz´ıvanou variantou diferenciáln´ı evoluce je DE/rand/1/bin. Jak uvid´ıme v následuj´ıc´ım zdrojovém textu, jej´ı naprogramován´ı v Matlabu je velmi snadné a lze ji zapsat na zhruba 25 ˇra´dc´ıch, pokud vyuˇzijeme funkce zrcad a nahvyb_expt popsané v pˇredcházej´ıc´ıch kapitolách a funkci rand1_bin, která vytvoˇr´ı nov´ y zkusm´ y bod y mutac´ı rand /1/ a binomick´ ym kˇr´ıˇzen´ım, zdrojov´ y text této funkce následuje. function y = rand1_bin(P, F, CR, expt) N = length(P(:,1)); d = length(P(1,:)) - 1; y = P(expt(1),1:d); vyb = nahvyb_expt(N, 3, expt); % three random points without expt r1 = P(vyb(1),1:d); r2 = P(vyb(2),1:d); r3 = P(vyb(3),1:d); u = r1 + F*(r2 - r3); change = find(rand(1,d) < CR); if isempty(change) % at least one element is changed change = 1 + fix(d*rand(1)); end y(change) = u(change);

50


function [x_star,fn_star,func_evals] = ... de_rand1_bin(fn_name, a, b, N, my_eps, max_evals, F, CR, shift) % input parameters: % a, b row vectors, limits of search space, a < b % N size of population % my_eps small positive value for stopping criterion % max_evals max. evals per one dimension of search space % F,CR mutation parameters % output: % fn_star the minimal function value found by MCRS % x_star global minimum found by search % func_evals number of func_evals for reaching stop d = length(a); P = zeros(N,d+1); for i=1:N P(i,1:d) = a + (b - a).*rand(1,d); P(i,d+1) = feval(fn_name, P(i,1:d)-shift); end % 0-th generation initialized fmax = max(P(:,d+1)); [fmin, indmin] = min(P(:,d+1)); func_evals = N; Q = P; while (fmax-fmin > my_eps) && (func_evals < d*max_evals) % main loop for i=1:N % in each generation y = rand1_bin(P ,F, CR, i); y = zrcad(y, a, b); % perturbation fy = feval(fn_name, y-shift); func_evals = func_evals + 1; if fy < P(i,d+1) % trial point y is good Q(i,:) = [y fy]; end end % end of generation P = Q; fmax = max(P(:,d+1)); [fmin,indmin] = min(P(:,d+1)); end % main loop - end x_star = P(indmin,1:d); fn_star = fmin;


51

Shrnut´ı: - Mutace jako pˇriˇcten´ı rozd´ılu vektor˚ u náhodnˇe vybran´ ych z populace - Kˇr´ıˇzen´ı v´ ymˇenou náhodnˇe vybran´ ych prvk˚ u vektoru - Binomické kˇr´ıˇzen´ı, exponenciáln´ı kˇr´ıˇzen´ı - V´ ybˇer lepˇs´ıho z dvojice (selekce turnajem)

Kontroln´ı ot´ azky: 1. Jaké jsou hlavn´ı odliˇsnosti ˇr´ızeného náhodného v´ ybˇeru a diferenciáln´ı evoluce? ˇ ım se liˇs´ı binomické a exponenciáln´ı kˇr´ıˇzen´ı? 2. C´ 3. Proˇc mus´ı b´ yt F ̸= 0? Co by se stalo, kdybychom zadali F < 0? 4. Co zp˚ usob´ı zadán´ı CR = 0 nebo CR = 1? Bude se pak liˇsit binomické kˇr´ıˇzen´ı od exponenciáln´ıho?


52

6.4

Experiment´ aln´ı porovn´ an´ı jednoduch´ ych evoluˇ cn´ıch algoritm˚ u

Pr˚ uvodce studiem: V této kapitole si ukáˇzeme experimentáln´ı porovnán´ı tˇr´ı jednoduch´ ych evoluˇcn´ıch algoritm˚ u na ˇsesti testovac´ıch funkc´ıch pˇri dvou u ´rovn´ıch dimenze u ´loh. Poˇc´ıtejte asi se dvˇema hodinami studia.

Máme porovnat u ´ˇcinnost tˇr´ı jednoduch´ ych algoritm˚ u na 6 testovac´ıch funkc´ıch z kapitoly 3.2 pro dimenze problému d = 2 a d = 10. U ˇctyˇrech funkc´ı s globáln´ım minimem v bodˇe x∗ = (0, 0, . . . , 0) uˇzijeme náhodn´ y posun popsan´ y v kapitole 3.2. Testovány jsou tyto algoritmy

ˇ ızené náhodné prohledáván´ı (CRS) se znáhodnˇenou reflex´ı simplexu podle • R´ (5.2), α = 4, velikost populace byla N = 10d. • Evoluˇcn´ı strategie ES(µ, λ) s adaptac´ı velikosti smˇerodatn´ ych odchylek podle pravidla jedné pˇetiny, µ = 10d a λ = 4 µ. • Diferenciáln´ı evoluce DE/rand/1/bin s parametry F = 0.5 a CR = 0.5. Velikost populace byla N = 5d.

Pro kaˇzdou u ´lohu bylo provedeno 100 opakován´ı. Ve vˇsech u ´lohách byla stejná podm´ınka ukonˇcen´ı, hledán´ı bylo ukonˇceno, kdyˇz rozd´ıl mezi nejvˇetˇs´ı a nejmenˇs´ı funkˇcn´ı hodnotou v populaci byl menˇs´ı neˇz 1e−6 nebo poˇcet vyhodnocen´ı funkc´ı dosáhl hodnotu 20000 d. Za u ´spˇeˇsné hledán´ı bylo povaˇzováno nalezen´ı bodu s funkˇcn´ı hodnotou menˇs´ı neˇz 1e − 4. Sledována byla spolehlivost algoritm˚ u v nacházen´ı správného ˇreˇsen´ı tj. poˇcet bˇeh˚ u, kdy nalezená minimáln´ı hodnota funkce je menˇs´ı neˇz 1e − 4 (v tabulkách je hodnota spolehlivosti uvedena ve sloupci R) a ˇcasová nároˇcnost algoritm˚ u vyjádˇrená jako pr˚ umˇern´ y poˇcet vyhodnocen´ı funkce potˇrebn´ y k dosaˇzen´ı podm´ınky ukonˇcen´ı, v tabulkách ve sloupci nfe. Tabulky následuj´ı dále v textu. Vu ´lohách s d = 2 mˇel nejvyˇsˇs´ı spolehlivost algoritmus CRS, nejrychleji podm´ınky ukonˇcen´ı dosahoval algoritmus DE/rand/1/bin se spolehlivost´ı srovnatelnou s algoritmem ES(µ, λ), ovˇsem tento algoritmus mˇel podstatnˇe vyˇsˇs´ı ˇcasové nároky neˇz diferenciáln´ı evoluce. Povˇsimnˇeme si v´ ysledk˚ u u Rosenbrockovy funkce. Tam s daleko nejmenˇs´ımi ˇcasov´ ymi nároky dosahoval spolehlivosti 100 (vˇsechny bˇehy nalezly správné ˇreˇsen´ı), zat´ımco u jin´ ych funkc´ı byly jeho ˇcasové nároky srovnatelné s ES(µ, λ) a v´ yraznˇe vyˇsˇs´ı oproti DE/rand/1/bin. Je to jedna z ukázek d˚ usledku No

6.4 Experimentáln´ı porovnán´ı jednoduch´ ych evoluˇcn´ıch algoritm˚ u

53

free lunch teorému, nˇejak´ y algoritmus je efektivn´ı jen pro nˇejakou skupinu problém˚ u, v jin´ ych problémech jsou jiné algoritmy efektivnˇejˇs´ı. d=2 Algoritmus →

CRS

Funkce↓

R

ES(µ, λ)

DE/rand/1/bin

R

R

nfe

87 2988 94

810

735 100 1715 98

421

nfe

ackley

100 2053

dejong1

100

nfe

griewank

80 3893

27 2284 70

1234

rastrig

94 1768

80 2158 84

654

947

83 8438 66

1554

schwefel0

94 1293

96 2360 76

636

pr˚ umˇer

95 1781

79 3324 81

885

rosen

100

d = 10 Algoritmus →

CRS

ES(µ, λ)

Funkce↓

R

nfe

R

ackley

19

12230

dejong1

99

griewank

19

15650

rastrig

1

rosen

75

schwefel0 pr˚ umˇer

R

nfe

6

32212 100

15756

7911 100

26068 100

7728

0

34072 100

37392

53392

0

30188 100

30920

53019

0 199700

0 200000 36

57034

nfe

DE/rand/1/bin

0

198601

0

34024

99

13644

18

59377

83

50673

Vu ´lohách s d = 10 je jasn´ ym v´ıtˇezem algoritmus DE/rand/1/bin, kter´ y v pˇeti testovac´ıch u ´lohách dosáhl témˇeˇr 100% spolehlivosti vˇetˇsinou i s nejmenˇs´ımi ˇcasov´ ymi nároky. Ovˇsem kompletnˇe selhal v hledán´ı minima Rosenbrockovy funkce, kde byl opˇet nej´ uspˇeˇsnˇejˇs´ı algoritmus CRS. Testovaná varianta evoluˇcn´ı strategie zcela selhala na 4 testovac´ıch funkc´ıch a dokonce i k ˇreˇsen´ı nejjednoduˇsˇs´ıho problému (De Jongova prvn´ı funkce) potˇrebovala nesrovnatelnˇe vˇetˇs´ı ˇcas neˇz ostatn´ı dva algoritmy. Experimenty ukázaly, ˇze i velmi jednoduché varianty stochastick´ ych algoritm˚ u jsou schopny vyˇreˇsit nˇekteré u ´lohy, ale rozhodnˇe nem˚ uˇzeme spoléhat na to, ˇze vyˇreˇs´ı


54

kaˇzdou u ´lohu. To ostatnˇe nem˚ uˇzeme oˇcekávat od ˇza´dného stochastického algoritmu. Témˇeˇr jistˇe by bylo moˇzno i pro jednoduché algoritmy v tomto testu naj´ıt jiné nastaven´ı hodnot ˇr´ıdic´ıch parametr˚ u, pro které by pak dosáhly lepˇs´ıch v´ ysledk˚ u. Ale takové hledán´ı vhodn´ ych hodnot ˇr´ıdic´ıch parametr˚ u pro ˇreˇsen´ y problém je ˇcasovˇe nároˇcné a jeho jedin´ y smysl je z´ıskat zkuˇsenosti v tomto nastavován´ı, které pak lze v omezené m´ıˇre vyuˇz´ıt i pro jiné u ´lohy. Vhodnˇejˇs´ı cesta je hledán´ı adaptivn´ıch algoritm˚ u, které si mˇen´ı hodnoty ˇr´ıdic´ıch parametr˚ u v pr˚ ubˇehu ˇreˇsen´ı konkrétn´ı u ´lohy samy, a tak dosahuj´ı vyˇsˇs´ı u ´ˇcinnosti pro ˇsirˇs´ı tˇr´ıdu optimalizaˇcn´ıch problém˚ u. Shrnut´ı: - Experimentáln´ı porovnán´ı u ´ˇcinnosti stochastick´ ych algoritm˚ u - Spolehlivost, ˇcasová nároˇcnost algoritmu - Popis experimentu, prezentace a interpretace v´ ysledk˚ u

Kontroln´ı ot´ azky: 1. Co je nutné specifikovat v popisu experimentu? 2. Které veliˇciny charakterizuj´ı u ´ˇcinnost algoritmu? 3. Lze dosaˇzené v´ ysledky nˇejak od˚ uvodnit?

55

7

Algoritmy modeluj´ıc´ı chov´ an´ı skupin

Mezi stochastick´ ymi algoritmy pro ˇreˇsen´ı problému globáln´ı optimalizace maj´ı v´ yznamné m´ısto algoritmy, které modeluj´ı sociáln´ı chován´ı skupiny (hejna, smeˇcky, kolonie) ˇziv´ ych jedinc˚ u. Nˇekdy jsou oznaˇcovány jako modely particle swarm intelligence [6], coˇz bychom mohli pˇreloˇzit jako kolektivn´ı inteligence skupiny. Tyto algoritmy se podobaj´ı evoluˇcn´ım algoritm˚ um v tom, ˇze jde o modely interakce v´ıce kandidát˚ u ˇreˇsen´ı, ale liˇs´ı se ve zp˚ usobu modelován´ı interakce mezi jedinci. V evoluˇcn´ıch algoritmech máme populaci, která se vyv´ıj´ı t´ım, ˇze prob´ıhá v´ ybˇer silnˇejˇs´ıch jedinc˚ u, jejich kˇr´ıˇzen´ı a mutace jejich genetického kódu. U algoritm˚ u popisovan´ ych v této kapitole máme skupinu (do jisté m´ıry je to analogie populace), ale jej´ı pohyb v prohledávaném prostoru ˇreˇsen´ı se ˇr´ıd´ı pravidly skupinového chován´ı. Mezi nejpopulárnˇejˇs´ı algoritmy této kategorie patˇr´ı Particle Swarm Optimization (PSO), Self-Organizing Migration Algorithm (SOMA) a algoritmus mravenˇc´ı kolonie (Ant Colony). S principy nˇekter´ ych takov´ ych algoritm˚ u se seznám´ıme v této kapitole a ukáˇzeme si také, jak mohou b´ yt implementovány v Matlabu.

7.1

Algoritmus PSO

Pr˚ uvodce studiem: V této kapitole se seznám´ıte s principy algoritmu PSO, kter´ y vznikl jako model chován´ı skupiny (hejna) ryb nebo pták˚ u. Poˇc´ıtejte asi se ˇctyˇrmi aˇz ˇsesti hodinami studia.

Algoritmus Particle Swarm Optimization (PSO), kter´ y jste se uˇz moˇzná poznali ´ v pˇredmˇetu Uvod do soft computingu [21], byl inspirován chován´ım hejn ryb a pták˚ u. Jedinci v takov´ ych hejnech maj´ı své individuáln´ı chován´ı a souˇcasnˇe jsou ovlivˇ nováni ostatn´ımi ˇcleny hejna. Pˇredstavme si, ˇze celé hejno tvoˇrené N jedinci hledá nˇejak´ y spoleˇcn´ y c´ıl a pohybuje se po etapách smˇerem k tomuto c´ıli. Zmˇena pozice jedince v dané etapˇe je dána jeho rychlost´ı (jak v´ıte ze základn´ı ˇskoly, rychlost je dráha uraˇzená za jednotku ˇcasu). Necht’ tedy etapa trvá jednu ˇcasovou jednotku a etapy postupného pohybu jedinc˚ u hejna oˇc´ıslujeme 0, 1, 2, . . .. Pozici i-tého jedince na zaˇca´tku následuj´ıc´ı etapy urˇc´ıme tak, ˇze k jeho souˇcasné pozici xi pˇriˇcteme dráhu, kterou v nynˇejˇs´ı etapˇe jedinec uraz´ı (tj. rychlost násobena délkou trván´ı etapy, ale uˇz jsme se dohodli, ˇze etapa trvá jednotkov´ y ˇcas, takˇze násoben´ı jedniˇckou m˚ uˇzeme vynechat). Pak nová pozice i-tého jedince v etapˇe t + 1 je dána vztahem xi (t + 1) = xi (t) + v i (t + 1).

(7.1)

56

´ Í SKUPIN 7 ALGORITMY MODELUJÍCÍ CHOVAN

Rychlost jedince v aktuáln´ı etapˇe je ovlivˇ nována jednak minulost´ı tohoto jedince, jednak minulost´ı celého hejna. Jak jedinec, tak i ostatn´ı ˇclenové hejna z´ıskali v minulosti nˇejaké znalosti o prohledávaném prostoru. Nejjednoduˇseji m˚ uˇzeme vyjádˇrit rychlost i-tého jedince v etapˇe t + 1 takto: v i (t + 1) = ω v i (t) + φ1 u1 ⊗ (pi − xi ) + φ2 u2 ⊗ (g − xi )

(7.2)

Koeficienty ω, φ1 , φ2 jsou vstupn´ı ˇr´ıdic´ı parametry algoritmu, u1 , u2 jsou vektory dimenze d, jejichˇz prvky jsou nezávislé náhodné hodnoty rovnomˇernˇe rozdˇelené na intervalu [0, 1), kter´ y v Matlabu m˚ uˇzeme vygenerovat pˇr´ıkazem rand(1,d). Symbol ⊗ znaˇc´ı násoben´ı dvou vektor˚ u po sloˇzkách, v Matlabu se pro tuto operaci uˇz´ıvá operátor .*. Vid´ıme, ˇze rychlost i-tého jedince v etapˇe t + 1 je tvoˇrena tˇremi sloˇzkami • Rychlost´ı i-tého jedince v minulé etapˇe v i (t), tedy setrvaˇcn´ a sloˇzka vyjadˇruj´ıc´ı smˇer, kter´ ym je jedinec rozbˇehnut. Tato rychlost je násobena váhou setrvaˇcˇ ım menˇs´ı je hodnota ω, t´ım v´ıce jedinec nosti ω, která se obvykle vol´ı ω < 1. C´ “zapom´ıná” svou rychlost z minulé etapy. • Druh´ y ˇclen φ1 u1 ⊗ (pi − xi ) je urˇcován minulost´ı jedince, pi je poloha bodu s nejmenˇs´ı funkˇcn´ı hodnotou, kterou i-t´ y jedinec naˇsel od zaˇca´tku prohledáván´ı. Hodnota pi b´ yvá v prac´ıch zab´ yvaj´ıc´ıch se PSO oznaˇcována jako personal best (pbest). Tento ˇclen odráˇz´ı to, co se jedinec ve své minulosti nauˇcil, oznaˇcuje se jako kognitivn´ı sloˇzka rychlosti, nˇekdy také jako nostalgie, protoˇze vyjadˇruje sklon jedince k návratu do nejlepˇs´ıho bodu své minulosti. • Tˇret´ı ˇclen φ2 u2 ⊗ (g − xi ) je urˇcován minulost´ı celého hejna, bod g (global best, gbest) je poloha bodu s nejmenˇs´ı funkˇcn´ı hodnotou, která byla nalezena hejnem od zaˇca´tku prohledáván´ı. Tento ˇclen je soci´ aln´ı sloˇzka rychlosti, nebot’ modeluje vliv spoleˇcnosti (hejna) na chován´ı jedince, tedy nˇeco jako podlehnut´ı vlivu okol´ı nebo davové psychóze. Geometrick´ y v´ yznam skládán´ı vektor˚ u podle vztah˚ u (7.1) a (7.2) je znázornˇen na následuj´ıc´ım obrázku. Povˇsimnˇeme si, ˇze vektor kognitivn´ı sloˇzky φ1 u1 ⊗ (pi − xi ) “nem´ıˇr´ı” pˇr´ımo do bodu pbest a podobnˇe ani vektor sociáln´ı sloˇzky φ2 u2 ⊗ (g − xi ) “nem´ıˇr´ı” pˇr´ımo do bodu gbest. Je to zp˚ usobeno t´ım, ˇze v obou tˇechto sloˇzkách je vektor rozd´ılu násoben náhodn´ ym vektorem, takˇze smˇer je ovlivˇ nován i náhodou a je pouze pˇribliˇznˇe do bod˚ u pbest, resp. gbest. Takto struˇcn´ y popis nám staˇc´ı k tomu, abychom mohli implementovat základn´ı nejjednoduˇsˇs´ı verzi algoritmu PSO. Pozice jedinc˚ u v nulté etapˇe se vygeneruj´ı náhodnˇe rovnomˇernˇe rozdˇelené v prohledávaném prostoru D. Podobnˇe i poˇca´teˇcn´ı rychlosti se vygeneruj´ı náhodnˇe. Zápis tohoto algoritmu v Matlabu následuje za obrázkem.

7.1 Algoritmus PSO

57

15

10

sociální slozka

gbest 5

v (t)

xi(t+1)

i

ω vi(t)

xi(t) kognitivní slozka

pbest

0

vi(t+1) −5 −5

0

5

10

15

function [xmin, fmin, func_evals, evals_near]=... pso1(fn_name, a, b, N, fi1, fi2, omega, ... my_eps, max_evals, shift, fnear) % d=length(a); P=zeros(N,d+1); % alokace pameti pro P v=zeros(N,d); % alokace pameti pro v vmax=(b-a)/6; % pro generovani pocatecnich rychlosti size(vmax) for i=1:N P(i,1:d)=a+(b-a).*rand(1,d); P(i,d+1)=feval(fn_name,P(i,1:d)-shift); v(i,:)= -vmax + 2*vmax.*rand(1,d); end % 0-th generation initialized func_evals=N; evals_near=0; pbest=P; % nejlepsi pozice jedincu na pocatku = pocatecni pozice [gbest_fun,indgbest]=min(pbest(:,d+1)); gbest=P(indgbest,:); fmax=max(P(:,d+1)); fmin=min(P(:,d+1)); while (fmax-fmin > my_eps) && (func_evals < d*max_evals) for i=1:N % in each generation

58


x=P(i,1:d); % pozice i-teho jedince v(i,:)=omega*v(i,:)+fi1*rand(1,d).*(pbest(i,1:d)-x)... +fi2*rand(1,d).*(gbest(1:d)-x); y = x + v(i,:); % nova pozice y=zrcad(y,a,b); % zabranuje opustit definicni obor D fy=feval(fn_name,y-shift); P(i,:)=[y,fy]; % zapis novou pozici a funkcni hodnotu func_evals=func_evals+1; if fy < pbest(i,d+1) % aktualizace pbest pbest(i,:)= [y,fy]; % prepis pbest end if fy < gbest(d+1) % aktualizace gbest gbest=[y, fy]; end end % end of generation fmax=max(P(:,d+1)); fmin=min(P(:,d+1)); if evals_near==0 && fmin
Tato základn´ı verze algoritmu PSO byla v pr˚ ubˇehu posledn´ıch let mnohokrát zdokonalována. Jedn´ım z takov´ ych krok˚ u bylo zaveden´ı koeficientu sevˇren´ı (coefficient of constriction). T´ım se drobnˇe zmˇen´ı vztah (7.2) pro urˇcen´ı rychlosti jedince v etapˇe t + 1: v i (t + 1) = χ [v i (t) + φ1 u1 ⊗ (pi − xi ) + φ2 u2 ⊗ (g − xi )] , (7.3) kde koeficient sevˇren´ı χ je urˇcen vztahem χ=

2κ √ , φ − 2 + φ2 − 4φ

kde φ = φ1 + φ2 a κ je vstupn´ı parametr ovlivˇ nuj´ıc´ı velikost koeficientu sevˇren´ı. Obvykle se vol´ı κ ≈ 1. Vzhledem k tomu, ˇze je nutné, aby φ > 4, je potˇreba tomuto poˇzadavku podˇr´ıdit volbu vstupn´ıch parametr˚ u φ1 a φ2 . Pro v literatuˇre doporuˇcované hodnoty φ1 = φ2 = 2.05 a κ = 1 vyjde hodnota koeficientu sevˇren´ı . χ = 0.73.

7.1 Algoritmus PSO

59

Dalˇs´ı modifikac´ı algoritmu PSO je jin´ y zp˚ usob v´ ypoˇctu sociáln´ı sloˇzky rychlosti jedince. Vytvoˇr´ıme model, kdy jedinec nen´ı ovlivˇ nován celou skupinou, ale jenom jej´ı ˇca´st´ı. Motivaci a neformáln´ı vysvˇetlen´ı k této u ´pravˇe algoritmu jste si mohli pˇreˇc´ıst v uˇcebn´ım textu [21], kdy rozpt´ ylená velká skupina turist˚ u hledá nejvyˇsˇs´ı vrchol Beskyd a jej´ı ˇclenové spolu komunikuj´ı vys´ılaˇckami omezeného dosahu. Informaci o v´ yˇsce a pozici jimi dosud nalezeného vrcholu mohou sdˇelit jen té ˇca´sti skupiny, která je v dosahu jejich vys´ılaˇcky. Zde si ukáˇzeme tzv. statickou topologii zvanou local best (lbest), kdy jedinci jsou rozdˇeleni do skupin, tyto skupiny se bˇehem celého prohledáván´ı nemˇen´ı a sociáln´ı sloˇzka rychlosti jedince nen´ı ovlivˇ nována bodem gbest, ale bodem lbest, coˇz je nejlepˇs´ı bod nalezen´ y skupinou, do které dan´ y jedinec patˇr´ı. Rychlost jedince v etapˇe t + 1 je pak urˇcena vztahem v i (t + 1) = χ [v i (t) + φ1 u1 ⊗ (pi − xi ) + φ2 u2 ⊗ (lj − xi )] ,

(7.4)

kde lj je lbest j-té skupiny, tj. souˇradnice nejlepˇs´ıho bodu dosud nalezeného skupinou j, do n´ıˇz patˇr´ı jedinec i. Ukáˇzeme si jednu z moˇznost´ı, jak m˚ uˇzeme implementovat statickou topologii lbest. Pro zadanou velikost hejna N a velikost skupin k urˇc´ıme poˇcet skupin s = ⌊N/k⌋, kde symbol ⌊·⌋ znaˇc´ı celou ˇca´st. Pokud N je dˇelitelné k beze zbytku, jsou vˇsechny skupiny stejnˇe velké, jinak posledn´ı skupina je vˇetˇs´ı. Na zaˇca´tku jsou jedinci do skupin pˇridˇeleni náhodnˇe (jelikoˇz souˇradnice jedinc˚ u hejna jsou inicializovány náhodnˇe, pak pˇriˇrazen´ı prvn´ıch k jedinc˚ u do prvn´ı skupiny atd. znamená náhodné rozdˇelen´ı skupin). Ze zkuˇsenosti vypl´ yvá, ˇze v pr˚ ubˇehu procesu hledán´ı se vˇetˇsinou skupiny shluknou tak, aby jedinci z téˇze skupiny si byli i geometricky bl´ızc´ı. Pokud bychom v tomto algoritmu zadali k = N (pˇr´ısnˇe vzato, staˇc´ı, aby k > N/2), pak by se topologie zmˇenila na gbest. Zdrojov´ y kód algoritmu PSO s koeficientem sevˇren´ı a se statickou topologi´ı lbest je uveden v následuj´ıc´ım v´ ypisu:

function [xmin, fmin, func_evals, evals_near]=... pso_lbest(fn_name, a, b, N, fi1, fi2, kappa, k,... my_eps, max_evals, shift, fnear) % % k je velikost skupin, pokud k=N, pak je to gbest % pocet_skupin=fix(N/k); % skupiny velikosti k, posledni skupina muze byt vetsi prislusnost=zeros(N,1); for ii=1:pocet_skupin-1 skup=ii+zeros(k,1);

60


prislusnost((ii-1)*k +1:(ii-1)*k+k) = skup; end prislusnost((pocet_skupin-1)*k +1:end) = pocet_skupin; d=length(a); P=zeros(N,d+1); v=zeros(N,d); vmax=(b-a)/6; % pro generovani pocatecnich rychlosti for i=1:N P(i,1:d)=a+(b-a).*rand(1,d); P(i,d+1)=feval(fn_name,P(i,1:d)-shift); v(i,:)= -vmax + 2*vmax.*rand(1,d); end % 0-th generation initialized pbest=P; % nejlepsi pozice jedincu na pocatku = pocatecni pozice lbest=zeros(pocet_skupin,d+1); for ii=1:pocet_skupin ind_skup= prislusnost==ii; pom=P(ind_skup,:); [tmp, ind_lmin]=min(pom(:,d+1)); lbest(ii,:)=pom(ind_lmin,:); end fmax=max(P(:,d+1)); fmin=min(P(:,d+1)); func_evals=N; evals_near=0; fi=fi1+fi2; chi=2*kappa / (fi-2+sqrt(fi*(fi-4))); % chi je koef. sevreni (constriction coeff.) while (fmax-fmin > my_eps) && (func_evals < d*max_evals) for i=1:N % in each generation x=P(i,1:d); % pozice i-teho jedince v(i,:)=chi*(v(i,:)+fi1*rand(1,d).*(pbest(i,1:d)-x)... +fi2*rand(1,d).*(lbest(prislusnost(i),1:d)-x)); y = x + v(i,:); % nova pozice y=zrcad(y,a,b); % zabranuje opustit definicni obor fy=feval(fn_name,y-shift); P(i,:)=[y,fy]; % zapis novou pozici a funkcni hodnotu func_evals=func_evals+1; if fy < pbest(i,d+1) % aktualizace pbest pbest(i,:)= [y,fy]; % prepis pozici a funkcni hodnotu pbest end

7.1 Algoritmus PSO

61

if fy < lbest(prislusnost(i),(d+1)) % aktualizace gbest lbest(prislusnost(i),:)=[y, fy]; end end % end of generation fmax=max(P(:,d+1)); fmin=min(P(:,d+1)); if evals_near==0 && fmin

62

PSO1 fce

nfe R

nfenear

PSO lbest std

nfe

R nfenear

15323 11160 32886 10

std

ackley

40000

7

dejong1

40000

9

griewank

40000

3

rastrig

40000

5

6700

2420 40000

rosen

40000

9

9584

8516 15548 10

2952 1325

schwefel0 40000

5

8396

3496 40000

2946 3098

1571

722 24776 10

19093 11498 40000

2260 1176 824

927

5

6716 8182

7

1777

7

788

Vid´ıme, ˇze algoritmus pso_lbest byl v hledán´ı globáln´ıho minima tˇechto funkc´ı u ´spˇeˇsnˇejˇs´ı. Zat´ımco pso1 ani jednou neukonˇcil bˇeh jinak neˇz dosaˇzen´ım maximáln´ıho zadaného poˇctu vyhodnocen´ı funkce (40000), pso_lbest ve tˇrech ze ˇsesti ˇreˇsen´ ych problém˚ u ukonˇcil u ´spˇeˇsnˇe hledán´ı pˇred dosaˇzen´ım maximáln´ıho zadaného poˇctu vyhodnocen´ı funkce. Nav´ıc, v tˇechto tˇrech problémech nalezl vˇzdy pˇrijatelnou aproximaci globáln´ıho minima, ve sloupci R je zaznamenám poˇcet takov´ ych u ´spˇeˇsn´ ych ˇreˇsen´ı. Rovnˇeˇz vid´ıme, ˇze u vˇsech testovan´ ych funkc´ı algoritmus pso_lbest byl v nalezen´ı dobrého ˇreˇsen´ı u ´spˇeˇsnˇejˇs´ı ˇcastˇeji neˇz pso1 (srovnej sloupce R) a pokud správné ˇreˇsen´ı nalezl, tak to bylo s menˇs´ı ˇcasovou nároˇcnost´ı (srovnej sloupce nfenear ), nav´ıc maj´ı tyto hodnoty menˇs´ı variabilitu neˇz u pso1, jak vid´ıme ze srovnán´ı hodnot jejich smˇerodatn´ ych odchylek ve sloupc´ıch oznaˇcen´ ych std. Z v´ ysledk˚ u tohoto srovnán´ı by mohl vzniknout dojem, ˇze i tak jednoduché verze algoritmu PSO jsou schopny ˇreˇsit obt´ıˇzné problémy globáln´ı optimalizace s vysokou u ´spˇeˇsnost´ı. Bohuˇzel tomu tak nen´ı. Vyzkouˇsejte si tyto algoritmy na stejn´ ych problémech, ale s vˇetˇs´ı dimenz´ı, tˇreba d = 5 nebo d = 10 a zjist´ıte, ˇze u ´spˇeˇsnost bude podstatnˇe niˇzˇs´ı. Podobnˇe nepˇr´ıznivˇe pro PSO dopadne i srovnán´ı s v´ ysledky jednoduch´ ych evoluˇcn´ıch algoritm˚ u v kapitole 6.4. Navzdory této skuteˇcnosti je algoritmus PSO spolu s diferenciáln´ı evoluc´ı povaˇzován za velmi nadˇejnou heuristiku pro ˇreˇsen´ı problém˚ u globáln´ı optimalizace a je v posledn´ıch létech intenzivnˇe studován. Byly navrˇzeny jeho mnohé sofistikované varianty, s nˇekter´ ymi se m˚ uˇzeme seznámit napˇr. ubˇeˇznˇe objevuj´ı v odborn´ ych ˇcasopisech. v monografii [6], dalˇs´ı se pr˚

7.1 Algoritmus PSO

63

Shrnut´ı: - Principy algoritmu PSO - Setrvaˇcná, kognitivn´ı a sociáln´ı sloˇzka rychlosti - Koeficient sevˇren´ı, topologie gbest a lbest

Kontroln´ı ot´ azky: 1. V ˇcem se liˇs´ı topologie gbest a lbest? Porovnejte v´ yhody a nev´ yhody tˇechto strategi´ı pˇri aplikac´ıch algoritmu PSO. 2. Porovnejte experimentálnˇe na ˇsesti testovac´ıch funkc´ıch z kap. 3.2 pro d = 2 a dvˇe r˚ uzné velikosti skupin v algoritmu pso lbest.


64

7.2

Algoritmus SOMA

Pr˚ uvodce studiem: V této kapitole se struˇcnˇe seznám´ıte s algoritmem SOMA, kter´ y modeluje chován´ı smeˇcky pronásleduj´ıc´ı koˇrist. Poˇc´ıtejte asi se tˇremi aˇz ˇctyˇrmi hodinami studia.

Algoritmus SOMA v roce 2000 navrhli Zelinka a Lampinen, takˇze jedn´ım z jeho autor˚ u je náˇs krajan. SOMA je zkratka pomˇernˇe dlouhého jména algoritmu (SelfOrganizing Migration Algorithm), které vˇsak vystihuje jeho základn´ı principy. Na algoritmus m˚ uˇzeme nahl´ıˇzet jako na jednoduch´ y model lov´ıc´ı smeˇcky. Ve smeˇcce je N jedinc˚ u a ti se pohybuj´ı (migruj´ı) po u ´zem´ı D tak, ˇze kaˇzd´ y jedinec v kaˇzdém migraˇcn´ım kole doskáˇce pˇr´ım´ ym smˇerem k v˚ udci smeˇcky, provˇeˇr´ı m´ısto kaˇzdého doskoku na této cestˇe a zapamatuje si nejlepˇs´ı j´ım nalezené m´ısto do dalˇs´ıho migraˇcn´ıho kola. V˚ udcem je vˇetˇsinou jedinec s nejlepˇs´ı hodnotou u ´ˇcelové funkce v celé smeˇcce a rovnˇeˇz kvalita m´ıst (bod˚ u v D) navˇst´ıven´ ych jedincem pˇri jeho skoc´ıch za v˚ udcem se posuzuje hodnotou u ´ˇcelové funkce. Pohyb migruj´ıc´ıho jedince nemus´ı prob´ıhat ve vˇsech dimenz´ıch prostoru D, ale jen v nˇekter´ ych dimenz´ıch, tedy v podprostoru s dimenz´ı menˇs´ı neˇz d. Migraci jedince k v˚ udci ukazuje následuj´ıc´ı obrázek. Jedinec z v´ ychoz´ıho bodu r 0 skáˇce skoky velikosti δ smˇerem k v˚ udci, pˇr´ıpadnˇe ho m˚ uˇze i o trochu pˇreskoˇcit. Pohyb jedince je bud’ ve smˇeru naznaˇceném plnou ˇsipkou, tedy v prostoru dimenze d nebo jen v prodprostoru menˇs´ı dimenze, tedy v nˇekterém ze smˇer˚ u vyznaˇcen´ ych teˇckovanˇe. Kromˇe vstupn´ıch parametr˚ u obvykl´ ych u vˇsech stochastick´ ych algoritm˚ u, tj. specifikace problému (funkce, prohledávan´ y prostor, podm´ınka ukonˇcen´ı) a velikosti smeˇcky (populace) N má algoritmus SOMA jeˇstˇe tˇri vstupn´ı parametry definuj´ıc´ı zp˚ usob pohybu jedince ve smˇeru za v˚ udcem.

7.2 Algoritmus SOMA

65

10 9

’

8

xL

7

x

L

6

(leader )

5 4

r 3

δ

2

r0

1 0

0

2

4

6

8

10

12

14

16

18

Pohyb jedince je ˇr´ızen tˇemito vstupn´ımi parametry (jejich oznaˇcen´ı ponecháme shodné s t´ım, které uˇz´ıvá ve své knize Zelinka [29] a jak se také ˇcasto oznaˇcuj´ı v jin´ ych publikac´ıch): • mass – relativn´ı velikost pˇreskoˇcen´ı v˚ udce ′

∥ xl − r 0 ∥ mass = , ∥ xl − r 0 ∥ kde r 0 je v´ ychoz´ı pozice jedince v daném migraˇcn´ım kole, xl je pozice v˚ udce, ′ ∥ xl − r 0 ∥ je norma vektoru (xl − r 0 ), podobnˇe ∥ xl − r 0 ∥ je norma vektoru ′ (xl − r 0 ). Tyto normy vektoru jsou vlastnˇe délky u ´seˇcek s koncov´ ymi body ′ r 0 , xl , resp. r 0 , xl . Geometrick´ y v´ yznam je také zˇrejm´ y z obrázku. Doporuˇceny jsou hodnoty mass ∈ [1.1; 3]. • step – urˇcuje velikost skok˚ u δ, tj. velikost smˇerového vektoru skoku podle následuj´ıc´ıho vztahu δ = step ∗ (xl − r 0 ), tedy ˇc´ım je hodnota parametru step menˇs´ı, t´ım podrobnˇeji je prostor na cestˇe k v˚ udci prozkoumáván. Doporuˇcené hodnoty jsou step ∈ [0.11; mass], autoˇri algoritmu SOMA rovnˇeˇz doporuˇcuj´ı volit hodnotu parametru step tak, aby 1/step nebylo celé ˇc´ıslo, tzn. aby jedinec na své cestˇe za v˚ udcem nenavˇst´ıvil zbyteˇcnˇe m´ısto obsazené a prozkoumané v˚ udcem. • prt – parametr pro urˇcen´ı smˇeru pohybu jedince za v˚ udcem, prt ∈ [0; 1]. Je to pravdˇepodobnost v´ ybˇeru dimenze prostoru D, ve které bude pohyb jedince prob´ıhat. Pro smˇer pohybu se vyberou náhodnˇe ty dimenze, pro které Ui <

66


prt, i = 1, 2, . . . , d. Hodnoty (U1 , U2 , . . . , Ud ) tvoˇr´ı náhodn´ y vektor velikosti d, jehoˇz sloˇzky jsou nezávislé náhodné veliˇciny rovnomˇernˇe rozdˇelené na intervalu [0, 1). Pˇri implementaci algoritmu se tyto hodnoty generuj´ı funkc´ı rand. Pokud by ˇza´dná sloˇzka Ui nesplˇ novala podm´ınku Ui < prt (napˇr. pˇri volbˇe prt = 0), vybere se jedna dimenze náhodnˇe z mnoˇziny {1, 2, . . . , d}. S podobn´ ym postupem jsme se setkali u binomického kˇr´ıˇzen´ı v diferenciáln´ı evoluci. Smˇer pohybu, tzn. dimenze, ve kter´ ych prob´ıhaj´ı skoky, se vol´ı vˇzdycky pro kaˇzdé migraˇcn´ı kolo jedince k v˚ udci znovu, zat´ımco velikost skok˚ u, tedy δ je stejná pro cel´ y bˇeh algoritmu SOMA. Dosud popsan´ y algoritmus je v literatuˇre oznaˇcován jako varianta SOMA all-toone, kdy v kaˇzdém migraˇcn´ım kole vˇsichni jedinci (kromˇe v˚ udce) putuj´ı za v˚ udcem. Implementace této varianty v Matlabu je uvedena v následuj´ıc´ım textu této kapitoly. Pohyb jedince smˇerem k v˚ udci je implementován funkc´ı go_ahead: % find the best point y with function value fy % on the way from r0 to leader function yfy = go_ahead(fn_name, r0, delta,... pocet_kroku, prt, a, b, shift) d = length(r0); tochange = find(rand(d,1) < prt); if isempty(tochange) tochange = 1 + fix(d*rand(1)); % at least one is changed end C = zeros(pocet_kroku,d+1); r = r0; for i=1:pocet_kroku r(tochange) = r(tochange) + delta(tochange); r = zrcad(r, a, b); fr = feval(fn_name, r - shift); C(i,:) = [r, fr]; end [fmin, indmin] = min(C(:,d+1)); yfy = C(indmin,:); % row with min. fy is returned

Tuto funkci pak volá algoritmus SOMA, kter´ y implementován jako funkce soma_to1, jej´ıˇz zdrojov´ y text následuje na dalˇs´ı stránce.

7.2 Algoritmus SOMA

% SOMA all to 1 function [x_star,fn_star,func_evals, evals_near] =... soma_to1(fn_name, a, b,... N,my_eps,max_evals,mass,step,prt,fnear, shift) % % fn_name function to be minized (M file) % a, b row vectors, limits of search space, a < b % N size of population % my_eps small positive value for stopping criterion % max_evals max. evals per one dimension of search space % mass,step,prt tunnig parameters of SOMA % d=length(a); P=zeros(N,d+1); for i=1:N P(i,1:d)=a+(b-a).*rand(1,d); P(i,d+1)= feval(fn_name,P(i,1:d) - shift); end % 0-th generation initialized P=sortrows(P,d+1); pocet_kroku=fix(mass/step); func_evals=N; evals_near=0; while ((P(N,d+1) - P(1,d+1)) > my_eps) ... && (func_evals < d*max_evals) % main loop leader=P(1,1:d); for i=2:N r0=P(i,1:d); % starting position of individual delta = (leader - r0)*step; % delta is vector (1,d) yfy = go_ahead(fn_name, r0, delta, pocet_kroku,... prt, a, b, shift); func_evals = func_evals + pocet_kroku; if yfy(d+1) < P(i,d+1) % trial point y is good P(i,:) = yfy; end end P=sortrows(P,d+1); % fmin is in the 1st row if evals_near == 0 && P(1,d+1) < fnear evals_near=func_evals; end end % main loop - end

67


68

x_star = P(1,1:d); fn_star = P(1,d+1); Dále jsou v literatuˇre zmiˇ novány dalˇs´ı dvˇe varianty algoritmu SOMA: • all-to-all – v kaˇzdém migraˇcn´ım kole se v˚ udcem stávaj´ı postupnˇe vˇsichni jedinci a ostatn´ı putuj´ı za nimi. Populace se aktualizuje novˇe nalezen´ ymi body aˇz po skonˇcen´ı migraˇcn´ıho kola, tzn. po dokonˇcen´ı putován´ı vˇsech ke vˇsem. Oproti pˇredchoz´ı variantˇe je prohledáván´ı prostoru d˚ ukladnˇejˇs´ı, ale ˇcasovˇe nároˇcnˇejˇs´ı, nebot’ nejsou tolik vyuˇzity znalosti z´ıskané z pˇredcházej´ıc´ıho hledán´ı o poloze bodu s nejmenˇs´ı funkˇcn´ı hodnotou ve smeˇcce. • all-to-all-adaptive – podobnˇe jako u pˇredchoz´ı varianty v kaˇzdém migraˇcn´ım kole se v˚ udcem stávaj´ı postupnˇe vˇsichni jedinci a ostatn´ı putuj´ı za nimi, ale jedinec v populaci je nahrazen bezprostˇrednˇe po ukonˇcen´ı jeho migrace za v˚ udcem, pokud pˇri této migraci byl nalezen bod s lepˇs´ı funkˇcn´ı hodnotou, takˇze následuj´ıc´ı jedinci v pr˚ ubˇehu tohoto migraˇcn´ıho kola uˇz putuj´ı k tomuto novému v˚ udci. Pˇ r´ıklad 7.1 Jako v pˇredchoz´ı kapitole, i zde uvedeme experimentáln´ı porovnán´ı algoritmu SOMA pˇri dvou r˚ uzn´ ych hodnotách ˇr´ıdic´ıho parametru prt na stejn´ ych ˇsesti testovac´ıch funkc´ıch uveden´ ych v kap. 3.2, dimenze problému opˇet byla d = 2. Ostatn´ı ˇr´ıdic´ı parametry algoritmu SOMA byly pro obˇe varianty shodné, mass = 2 a step = 0.11. Pro kaˇzdou variantu a problém bylo provedeno 10 opakován´ı. V obou variantách byla velikost smeˇcky zvolena N = 20, podm´ınka ukonˇcen´ı definovaná vstupn´ımi parametry max_evals= 20000 a my_eps= 1e − 4. Podobnˇe jako v pˇredchoz´ı kapitole byla sledována nejen celková ˇcasová nároˇcnost daná poˇctem vyhodnocen´ı funkce nfe potˇrebn´ ych k dosaˇzen´ı podm´ınky ukonˇcen´ı, ale i poˇcet vyhodnocen´ı funkce potˇrebné k tomu, aby se nalezené minimum funkce pˇribl´ıˇzilo hodnotˇe v globáln´ım minimu testované funkce tak, ˇze se liˇs´ı o ménˇe neˇz 1 × 10−4 . Tento poˇcet vyhodnocen´ı je v tabulce oznaˇcen jako nfenear a je to pr˚ umˇerná hodnota ˇcasové nároˇcnosti u ´spˇeˇsn´ ych ˇreˇsen´ı, smˇerodatné odchylky jsou ve sloupc´ıch oznaˇcn´ ych std. V´ ysledky porovnán´ı jsou v následuj´ıc´ı tabulce, poˇcty vyhodnocen´ı funkce nfe jsou pr˚ umˇery z 10 opakován´ı.

7.2 Algoritmus SOMA

69

prt = 0.5 fce

nfe

R

nfenear

prt = 1 std

nfe

R

nfenear

std

ackley

3987 10

2380 340

5697

2

2414

967

dejong1

2756 10

1183 177

1525 10

1217

242

griewank

7202 10

2619 825

6689

1

2072

rastrig

3474 10

2004 270 12366

3

1730

rosen

6586 10

2072 773

4740

5

1798 1065

1987 242

5253

2

5663 1209

schwefel0 3508

8

592

Vid´ıme, ˇze varianta s prt = 0.5 byla v hledán´ı globáln´ıho minima tˇechto funkc´ı spolehlivˇejˇs´ı (srovnej sloupce R) a nav´ıc, v polovinˇe testovan´ ych funkc´ı i rychlejˇs´ı.

Shrnut´ı: - Model pohybu jedince za v˚ udcem a parametry mass, step, prt - Varianta all-to-one - Varianta all-to-all - Varianta all-to-all-adaptive

Kontroln´ı ot´ azky: 1. Proˇc má b´ yt parametr step volen tak, aby 1/step nebylo celé ˇc´ıslo? 2. Jak se liˇs´ı varianta all-to-all od all-to-all-adaptive?

ˇ ZE ˇ STRATEGIÍ 8 ADAPTACE POMOCÍ SOUTE

70

8

Adaptace pomoc´ı soutˇ eˇ ze strategi´ı

Pr˚ uvodce studiem: Tato kapitola je vˇenována jedné z moˇznost´ı, jak implementovat adaptivn´ı stochastické algoritmy, tj. algoritmy, u kter´ ych nen´ı nutné pracnˇe nastavovat jejich ˇr´ıdic´ı parametry podle právˇe ˇreˇseného problému. Poˇc´ıtejte asi s pˇeti hodinami studia.

V minul´ ych kapitolách jsme se seznámili s nˇekolika stochastick´ ymi algoritmy a v´ıme uˇz, ˇze v r˚ uzn´ ych algoritmech se uˇz´ıvá r˚ uzná strategie hledán´ı nebo v jednom algoritmu uˇz´ıvaj´ıc´ım jednu strategii hledán´ı m˚ uˇze m´ıt tato strategie r˚ uznˇe zvolené hodnoty jej´ıch ˇr´ıdic´ıch parametr˚ u. Pokud v jednom stochastickém algoritmu nab´ıdneme nˇekolik takov´ ych stˇr´ıdaj´ıc´ıch se strategi´ı a mezi nimi vyb´ıráme podle jejich dosavadn´ı u ´spˇeˇsnosti hledán´ı s t´ım, ˇze preferujeme ty u ´spˇeˇsnˇejˇs´ı strategie, dáváme tak algoritmu moˇznost pˇrednostnˇe vyb´ırat strategii vhodnou pro právˇe ˇreˇsen´ y problém. Podrobnˇeji vysvˇetl´ıme principy této soutˇeˇze na zobecnˇeném algoritmu ˇr´ızeného náhodného prohledáván´ı (CRS). generuj populaci P , tj. N bod˚ u náhodnˇe v D; repeat najdi xmax ∈ P , f (xmax ) ≥ f (x), x ∈ P ; repeat uˇzij lokáln´ı heuristiku k vygenerov´ an´ı nového bodu y ∈ D; until f (y) < f (xmax ); xmax := y; until podm´ınka ukonˇcen´ı; Nov´ y bod y se v klasické variantˇe algoritmu CRS generuje reflex´ı simplexu, viz kapitola (5). Samozˇrejmˇe je vˇsak moˇzné ke generován´ı bodu y uˇz´ıt i jinou lokáln´ı heuristiku, jak uˇz bylo zm´ınˇeno v kapitole o algoritmu CRS, napˇr. znáhodnˇenou reflexi podle (5.2). Bod y vˇsak lze generovat jakoukoli jinou heuristikou. Volbou lokáln´ı heuristiky vol´ıme strategii hledán´ı v algoritmu CRS. Odtud uˇz je jen kr˚ uˇcek k uˇzit´ı v´ıce lokáln´ıch heuristik v rámci jednoho algoritmu CRS a tyto heuristiky stˇr´ıdat. Pˇri tomto stˇr´ıdán´ı budeme preferovat ty heuristiky, které byly v dosavadn´ım pr˚ ubˇehu hledán´ı u ´spˇeˇsnˇejˇs´ı. V takovém algoritmu heuristiky soutˇeˇz´ı o to, aby byly vybrány pro generován´ı nového bodu. Protoˇze dáváme pˇrednost u ´spˇeˇsnˇejˇs´ım heuristikám, je nadˇeje, ˇze bude vybrána vhodná strategie pro ˇreˇsen´ y problém a algoritmus bude konvergovat rychleji. Takto vyuˇz´ıváme uˇcen´ı algo-

71

ritmu v pr˚ ubˇehu ˇreˇsen´ı daného problému. Vyb´ıráme pˇrednostnˇe tu heuristiku, která má vˇetˇs´ı fitness ohodnocovanou dosavadn´ı u ´spˇeˇsnost´ı. Mˇejme tedy v algoritmu k dispozici h strategi´ı (v pˇr´ıpadˇe algoritmu CRS lokáln´ıch heuristik) a v kaˇzdém kroku ke generován´ı nového bodu vyb´ıráme náhodnˇe i-tou heuristiku s pravdˇepodobnost´ı qi , i = 1, 2, . . . , h. Na zaˇcátku vˇsechny pravdˇepodobnosti nastav´ıme na shodné hodnoty qi = 1/h. Pravdˇepodobnosti qi mˇen´ıme v závislosti na u ´spˇeˇsnosti i-té heuristiky v pr˚ ubˇehu vyhledávac´ıho procesu. Heuristiku povaˇzujeme za u ´spˇeˇsnou, kdyˇz generuje takov´ y nov´ y bod y, ˇze f (y) < f (xmax )), tj. nov´ y bod y je zaˇrazen do populace. Pravdˇepodobnosti qi m˚ uˇzeme vypoˇc´ıtat jako relativn´ı ˇcetnost u ´spˇech˚ u i-té heuristiky. Pokud ni je dosavadn´ı poˇcet u ´spˇech˚ u i-té heuristiky, pravdˇepodobnost qi je u ´mˇerná tomuto poˇctu u ´spˇech˚ u qi = ∑h

ni + n0

j=1 (nj

+ n0 )

,

(8.1)

kde n0 > 0 je vstupn´ı parametr algoritmu. Nastaven´ım n0 > 1 zabezpeˇc´ıme, ˇze jeden náhodn´ yu ´spˇech heuristiky nevyvolá pˇr´ıliˇs velkou zmˇenu v hodnotˇe qi . Algoritmus uˇz´ıvaj´ıc´ı k hodnocen´ı u ´spˇeˇsnosti heuristik vztah (8.1) je jednou z moˇzn´ ych variant hodnocen´ı fitness heuristiky. Jinou moˇznost´ı, jak ohodnotit u ´spˇeˇsnost heuristiky, je váˇzit u ´spˇeˇsnost relativn´ı zmˇenou v hodnotˇe funkce. Váha wi se urˇc´ı jako wi =

fmax − max(f (y), fmin ) , fmax − fmin

(8.2)

kde fmax a fmin jsou nejvˇetˇs´ı, resp. nejmenˇs´ı funkˇcn´ı hodnota v populaci. Hodnoty wi jsou tedy v intervalu (0, 1) a pravdˇepodobnost qi se pak vyhodnot´ı jako qi = ∑h

Wi + w0

j=1

(Wj + w0 )

,

(8.3)

kde Wi je souˇcet vah wi v pˇredcházej´ıc´ım hledán´ı a w0 > 0 je vstupn´ı parametr algoritmu. Aby se zabránilo degeneraci rozdˇelen´ı pravdˇepodobnost´ı qi napˇr. pˇr´ıliˇsnou preferenc´ı heuristiky, které byla u ´spˇeˇsná na zaˇca´tku hledán´ı, lze zavést vstupn´ı parametr δ a klesne-li kterákoli hodnota qi pod hodnotu δ, jsou pravdˇepodobnosti v´ ybˇeru heuristik nastaveny na jejich poˇcáteˇcn´ı hodnoty qi = 1/h. Dalˇs´ı prostor pro tvorbu r˚ uzn´ ych adaptivn´ıch variant algoritmu je volba soutˇeˇz´ıc´ıch strategi´ı (v pˇr´ıpadˇe algoritmu CRS lokáln´ıch heuristik), které jsou zaˇrazeny mezi h soutˇeˇz´ıc´ıch lokáln´ıch heuristik. Pˇr´ıklady moˇzn´ ych lokáln´ıch heuristik jsou reflexe a

72


znáhodnˇená reflexe v simplexu popsané uˇz v kapitole 5. Dalˇs´ı moˇznost´ı jsou lokáln´ı heuristiky vycházej´ıc´ı z diferenciáln´ı evoluce, kdy bod u se generuje jako v difereciáln´ı evoluci a nov´ y vektor y vznikne kˇr´ıˇzen´ım vektoru u a vektoru x n´ ahodnˇe vybraného z populace. Heuristiky inspirované evoluˇcn´ı strategi´ı mohou genererovat nov´ y bod y podle následuj´ıc´ıho pravidla yi = xi + U, i = 1, 2, . . . , d, kde xi je i-tá souˇradnice náhodnˇe vybraného bodu x z populace a U je náhodná veliˇcina, U ∼ N (0, σi2 ). Zat´ımco v evoluˇcn´ı strategii je hodnota parametru σi2 adaptována v pr˚ ubˇehu procesu vyhledáván´ı podle tzv. pravidla 1/5 u ´spˇechu, zde m˚ uˇze b´ yt aktuáln´ı hodnota tohoto parametru odvozována z aktuáln´ı populace, napˇr. jako σi = k (max xi − min xi ) + ε, x∈P x∈P

i = 1, 2, . . . , d,

operátory max, min znamenaj´ı nejvˇetˇs´ı, resp. nejmenˇs´ı hodnotu i-té souˇradnice v aktuáln´ı populaci P , vstupn´ı parametr ε > 0 (v implementaci m˚ uˇzeme uˇz´ıt napˇr. −4 ε = 1 × 10 ) zabezpeˇcuje, ˇze hodnota σi je kladná, k je vstupn´ı ˇr´ıdic´ı parametr heuristiky, k > 0. Lokáln´ı heuristikou v adaptivn´ım algoritmu CRS m˚ uˇze b´ yt samozˇrejmˇe také náhodné generován´ı bodu y ze spojitého rovnomˇerného rozdˇelen´ı na D, které se uˇz´ıvá ve slepém náhodném prohledáván´ı. T´ım m˚ uˇzeme dokonce dosáhnout toho, ˇze u takového algoritmu lze teoreticky dokázat, ˇze je konvergentn´ı (o teoretické konvergenci algoritmu viz kapitola 4). Algoritmus se soutˇeˇz´ıc´ımi strategiemi je konvergentn´ı, kdyˇz (a) mezi soutˇeˇz´ıc´ımi strategiemi je alespoˇ n jedna zaruˇcuj´ıc´ı kladnou pravdˇepodobnost vygenerován´ı nového bodu y ∈ S, kde S je libovolná otevˇrená podmnoˇzina D maj´ıc´ı Lebesgueovu m´ıru λ(S) > 0. Takovou strategi´ı v pˇr´ıpadˇe CRS se soutˇeˇz´ıc´ımi lokáln´ımi heuristikami je náhodné generován´ı bodu y ze spojitého rovnomˇerného rozdˇelen´ı na D uˇz´ıvané ve slepém náhodném prohledáván´ı. (b) hodnoty pravdˇepodobnosti v´ ybˇeru takové strategie qk v kaˇzdém z k krok˚ u hle∑∞ dán´ı tvoˇr´ı divergentn´ı ˇradu ( k=1 qk = ∞). To je zajiˇstˇeno volbou vstupn´ıho parametru δ > 0. (c) nejlepˇs´ı bod staré populace vˇzdy pˇrecház´ı do nové populace (v algoritmu se uˇz´ıvá tzv. elitismus, coˇz v CRS je dodrˇzováno). Moˇzná, ˇze ve vás text této kapitoly vzbudil dojem, ˇze implementace algoritm˚ u se soutˇeˇz´ı strategi´ı je obzvláˇst’ obt´ıˇzná. Tento dojem je ale zbyteˇcn´ y. V´ ybˇer strategie podle jej´ı dosavadn´ı u ´spˇeˇsnosti lze snadno realizovat funkc´ı ruleta, kterou známe uˇz

73

z kapitoly o genetick´ ych algoritmech. Tato funkce dostane na vstupu vektor ˇcetnost´ı dosavadn´ı u ´spˇeˇsnosti jednotliv´ ych strategi´ı a vrát´ı ˇc´ıslo “vylosované” strategie spolu s minimáln´ı hodnotou pravdˇepodobnosti jednotliv´ ych strategi´ı. Implementace této funkce je v následuj´ıc´ım v´ ypisu zdrojového textu. function [hh, p_min]=ruleta(cutpoints) % returns an integer from [1, length(cutpoints)] % with probability proportional % to cutpoints(i)/ summa cutpoints h = length(cutpoints); ss = sum(cutpoints); p_min = min(cutpoints)/ss; cp = zeros(1, h); cp(1) = cutpoints(1); for i = 2:h cp(i) = cp(i-1) + cutpoints(i); end cp = cp/ss; hh = 1 + fix(sum(cp my_eps)&(func_evals < d*max_evals) %main loop [hh p_min]=ruleta(ni); if p_min<delta ni=zeros(1,h)+n0; % reset end switch hh case 1 y = lokalni heuristika1 case 2 y = lokalni heuristika2 ..... case 8 y = lokalni heuristika8 end y=zrcad(y,a,b); % perturbation fy=feval(fn_name,y); func_evals=func_evals+1;

74


if fy < fmax % trial point y is good ni(hh)=ni(hh)+1; % zmena prsti qi ..... end end % main loop - end

Podobnˇe lze snadno navrhnout a implementovat adaptivn´ı verze diferenciáln´ı evoluce, kde bude soutˇeˇzit nˇekolik r˚ uzn´ ych strategi´ı. Strategie se mohou liˇsit druhem mutace, typem kˇr´ıˇzen´ı nebo r˚ uzn´ ym nastaven´ım ˇr´ıdic´ıch parametr˚ u F a CR. R˚ uzné adaptivn´ı verze diferenciáln´ı evoluce byly v posledn´ıch létech intenzivnˇe zkoumány a byly navrˇzeny verze, které se osvˇedˇcily jak na testovac´ıch problémech, tak i v aplikac´ıch. D˚ uvodem pro hledán´ı adaptivn´ıch algoritm˚ u je pˇredevˇs´ım to, aby se pˇri jejich uˇzit´ı usnadnilo nastavován´ı ˇr´ıdic´ıch parametr˚ u, pˇr´ıpadnˇe vylouˇcilo u ´plnˇe, a pˇritom aby takové verze algoritm˚ u umoˇzn ˇovaly ˇreˇsit ˇsirokou skupinu problém˚ u s dobrou spolehlivost´ı a v pˇrijatelném ˇcase. Samozˇrejmˇe, ˇze existuj´ı i jiné zp˚ usoby adaptace ˇr´ıdic´ıch parametr˚ u, neˇz je uveden´ y mechanismus soutˇeˇze strategi´ı. Pˇripomeneme si algoritmus diferenciáln´ı evoluce: generuj poˇcáteˇcn´ı populaci P (N bod˚ u náhodnˇe v D) repeat for i := 1 to N do vytvoˇr vektor y mutac´ı a kˇr´ıˇzen´ım if f (y) < f (xi ) then do Q zaˇrad’ y else do Q zaˇrad’ xi endfor P := Q until podm´ınka ukonˇcen´ı Chceme-li implementovat diferenciáln´ı evoluci se soutˇeˇz´ı r˚ uzn´ ych strategi´ı, pak jenom potˇrebujeme zaˇr´ıdit, aby se v pˇr´ıkazu “vytvoˇr vektor y mutac´ı a kˇr´ıˇzen´ım” vyb´ıralo ruletou z h nab´ıdnut´ ych strategi´ı liˇs´ıc´ıch se typem mutace nebo kˇr´ıˇzen´ı nebo hodnotami parametr˚ u F a CR. Implementovat to m˚ uˇzeme pomoc´ı pˇr´ıkazu switch a mus´ıme doplnit naˇc´ıtán´ı hodnot u ´spˇeˇsnosti jednotliv´ ych strategi´ı pro v´ ypoˇcet hodnot pravdˇepodobnosti qi analogicky, jako bylo ukázáno u algoritmu CRS. Pˇ r´ıklad 8.1 V´ yznam a u ´ˇcinnost adaptivn´ıch algoritm˚ u ilustruj´ı v´ ysledky experimentáln´ıho srovnán´ı tˇr´ı takov´ ych algoritm˚ u na 6 testovac´ıch funkc´ıch z kapitoly 3.2. Dimenze problému byla pro vˇsechny funkce d = 10 a pro kaˇzdou funkci bylo provedeno 100 opakován´ı. Ve vˇsech u ´lohách byla stejná podm´ınka ukonˇcen´ı, hledán´ı bylo ukonˇceno, kdyˇz rozd´ıl mezi nejvˇetˇs´ı a nejmenˇs´ı funkˇcn´ı hodnotou v populaci byl

75

menˇs´ı neˇz 1e−6 nebo poˇcet vyhodnocen´ı funkce dosáhl hodnotu 20000d. Za u ´spˇeˇsné hledán´ı bylo povaˇzováno nalezen´ı bodu s funkˇcn´ı hodnotou menˇs´ı neˇz 1e − 4. Porovnávány byly tyto algoritmy: • jde [5] – adaptivn´ı diferenciáln´ı evoluce DE/rand/1/bin, ve které se hodnoty parametr˚ u F a CR nastavuj´ı evoluc´ı v pr˚ ubˇehu procesu hledán´ı. • b6e6rl [26] – diferenciáln´ı evoluce, ve které se strategie hledán´ı adaptuj´ı soutˇeˇz´ı. Jako mutace je ve vˇsech strategi´ıch uˇzita /randrl/1/, soutˇeˇz´ı 6 r˚ uzn´ ych dvojic nastaven´ı hodnot parametr˚ u F a CR pro binomické kˇr´ıˇzen´ı a 6 r˚ uzn´ ych dvojic nastaven´ı hodnot parametr˚ u F a CR pro exponenciáln´ı kˇr´ıˇzen´ı. • crs4hc [25] – CRS se soutˇeˇz´ı 4 lokáln´ıch heuristik. Jedin´ y ˇr´ıdic´ı parametr, kter´ y je nutno zadávat, je velikost populace. Pro obˇe varianty diferenciáln´ı evoluce bylo NP = 40, pro crs4hc NP = 100. Parametry ˇr´ıd´ıc´ı adaptaci byly ve vˇsech algoritmech ponechány na jejich implicitn´ım nastaven´ı doporuˇcovaném autory.


76

V´ ysledky jsou v následuj´ıc´ı tabulce:

Algoritmus →

jde

Funkce↓

R

b6e6rl nfe

R

crs4hc

nfe

R

nfe

ackley

100 16316

100 14899

100

18207

dejong1

100

100

7109

100

10017

98 31897

98 26099

75

69602

rastrig

100 18664

100 14058

rosen

100 68774

99 24081

100

18532

99 13296

100 11770

100

33685

99.5 26219 99.5 16336 83.3

56066

griewank


8367

25 186353

Porovnejte tyto v´ ysledky s v´ ysledky, které byly dosaˇzeny jednoduch´ ymi neadaptivn´ımi evoluˇcn´ımi algoritmy pro stejné problémy, které jsou uvedeny v kapitole 6.4. Pro pohodlnˇejˇs´ı porovnán´ı je zde tabulka s v´ ysledky zopakována. Vid´ıme, ˇze adaptivn´ı algoritmy ve vˇetˇsinˇe problém˚ u nacházej´ı dobré pˇribl´ıˇzen´ı globáln´ımu minimu s podstatnˇe vyˇsˇs´ı spolehlivost´ı a vˇetˇsinou i v kratˇs´ım ˇcase.

d = 10 Algoritmus →

CRS

ES(µ, λ)

Funkce↓

R

nfe

R

ackley

19

12230

dejong1

99

griewank

19

15650

rastrig

1

rosen

75


R

nfe

6

32212 100

15756

7911 100

26068 100

7728

0

34072 100

37392

53392

0

30188 100

30920

53019

0 199700

0 200000 36

57034

nfe

DE/rand/1/bin

0

198601

0

34024

99

13644

18

59377

83

50673

77

Shrnut´ı: - Adaptace algoritmu soutˇeˇz´ı strategi´ı hledán´ı - Soutˇeˇz lokáln´ıch heuristik v algoritmu CRS, pravidla pro posuzován´ı dosavadn´ı u ´spˇeˇsnosti - Podm´ınky konvergence algoritmu - Soutˇeˇz strategi´ı v diferenciáln´ı evoluci

Kontroln´ı ot´ azky: 1. Co rozum´ıme pojmem “lokáln´ı heuristika” v algoritmu CRS? 2. V ˇcem se liˇs´ı pravdˇepodobnosti v´ ybˇeru strategie podle vztahu (8.1) a (8.3)? 3. Jakou kombinaci strategi´ı navrhnete, aby algoritmus CRS se soutˇeˇz´ı lokáln´ıch heuristik splˇ noval podm´ınku konvergence? ˇ ım se mohou liˇsit soutˇeˇz´ıc´ı strategie v adaptivn´ım algoritmu diferenciáln´ı 4. C´ evoluce?

78

9

Literatura

[1] Ali M. M., Törn A., Population set based global optimization algorithms: Some modifications and numerical studies, Computers and Operations Research, 31, 2004, 1703–1725. [2] Bäck, T., Evolutionary Algorithms in Theory and Practice. Oxford University Press, New York, 1996. [3] Bäck, T., Schwefel, H.P., An Overview of Evolutionary Algorithms for Parameter Optimization, Evolutionary Computation 1, 1993, 1–23. [4] Bäck, T., Fogel, D.B., Michalewicz, Z. (eds). Evolutionary Computation 2 Advanced Algorithms and Operators, Institute of Physics Publishing, Bristol and Philadelphia, 2000. ˇ [5] Brest, J., Greimer, S., Boˇskoviˇc, B., Mernik, M., Zumer, V., Self-adapting Control Parameters in Differential Evolution: A Comparative Study on Numerical Benchmark Problems. IEEE Transactions on Evolutionary Computation, 10, 2006, 646–657. [6] Engelbrecht A. P., Computatinal Intelligence: An Introduction. Chichester: John Wiley & sons, 2007. [7] Feoktistov V., Differential Evolution in Search of Sotution. Springer, 2006. [8] Goldberg, D. E., Genetic Algorithms in Search, Optimization, and Machine Learning. Reading, Addison Wesley, 1989. [9] Holland J. H., Adaptation in Natural and Artificial Systems, The University of Michigan Press, Ann Arbor, MI, 1975. [10] Hynek J., Genetické algoritmy a genetické programov´ an´ı, Grada, 2008 [11] Kaelo P., Ali M. M., A numerical study of some modified differential evolution algorithms, European J. Operational Research, 169, 2006, 1176–1184. [12] Kvasniˇcka, V., Posp´ıchal, J., Tiˇ no, P., Evoluˇcné algoritmy. STU, Bratislava, 2000. [13] Michalewicz, Z., Genetic Algorithms + Data Structures = Evolution Programs. Berlin, Springer Verlag, 1992. [14] Michalewicz, Z., Fogel, D.B. How to Solve It: Modern Heuristics. SpringerVerlag, Berlin Heidelberg New York, 2000. ˇ Plzeˇ [15] M´ıka, S., Matematická optimalizace. vydavatelstv´ı ZCU, n, 1997 [16] Nelder, J.A., Mead, R., A Simplex Method for Function Minimization, Computer J., 7(1), 1964, 308–313. [17] Price, K. V., Storn, R. and Lampinen, J., Differential Evolution: A Practical Approach to Global Optimization. Springer, 2005.

79

[18] Price, W. L., A Controlled Random Search Procedure for Global Optimization. Computer J., 20, 1977, 367–370. [19] Spall J. C., Introduction to Stochastic Search and Optimization, WileyIntersience, 2003. [20] Storn, R., Price, K., Differential Evolution – a Simple and Efficient Heuristic for Global Optimization. J. Global Optimization, 11, 1997, 341–359. ˇ epniˇcka M., Vavˇr´ıˇcková L., Uvod ´ [21] Stˇ do soft computingu, Uˇcebn´ı text pro distanˇcn´ı studium, Ostravská universita, 2010 ˇ [22] Törn, A., Zilinskas A., Global Optimization, Lecture Notes in Computer Science, No. 350, Springer, 1989. [23] Tvrd´ık, J., Evoluˇcn´ı algoritmy. uˇcebn´ı text pro kombinované studium, Ostravská universita, 2010. http://albert.osu.cz/tvrdik/down/vyuka.html [24] Tvrd´ık, J., Pavliska V., Bujok P., Z´ aklady modelov´ an´ı v MATLABU. uˇcebn´ı text pro distanˇcn´ı a kombinované studium, Ostravská universita, 2010. [25] Tvrd´ık J., Kˇriv´ y I., Miˇs´ık L., Adaptive Population-based Search: Application to Estimation of Nonlinear Regression Parameters, Computational Statistics and Data Analysis 52, 2007, 713-724. [26] Tvrd´ık J., Self-adaptive Variants of Differential Evolution with Exponential Crossover. Analele Universitatii de Vest, Timisoara. Seria MatematicaInformatica, 47, 2009, 151–168. Reprint available online: http://albert.osu.cz/tvrdik/down/global_optimization.html [27] Wolpert, D. H., Macready, W.G., No Free Lunch Theorems for Optimization, IEEE Transactions on Evolutionary Computation 1 (1), 1997, 67–82. [28] Zaharie D., Influence of crossover on the behavior of Differential Evolution Algorithms. Applied Soft Computing, 9, 2009, 1126–1138. [29] Zelinka, I., Umˇel´ a inteligence v problémech globáln´ı optimalizace. BEN, Praha, 2002 ˇ [30] Zelinka I., Oplatková Z. Seda M., Oˇsmera P. Vˇcelaˇr F., Evoluˇcn´ı výpoˇcetn´ı techniky, BEN, 2009

80

WWW stránky: http://albert.osu.cz/oukip/optimization/ Knihovna algoritm˚ u v Matlabu a C++ k volnému vyuˇzit´ı. http://www.mat.univie.ac.at/~neum/glopt.html Obsáhlé stránky o globáln´ı optimalizaci a evoluˇcn´ıch algoritmech, mnoho dalˇs´ıch uˇziteˇcn´ ych odkaz˚ u na souvisej´ıc´ı témata vˇcetnˇe matematiky a statistiky, adresy stránek mnoha autor˚ u ˇclánk˚ u a knih s touto problematikou atd. http://www.cs.sandia.gov/opt/survey/main.html pˇrehled, odkazy na základn´ı ˇclánky z 90. let a dˇr´ıvˇejˇs´ı doby, od roku 1997 nejsou tyto stránky aktualizovány

Stochastické algoritmy pro globální optimalizaci

Recommend Documents