UČEBNÍ TEXTY OSTRAVSKÉ UNIVERZITY. Přírodovědecká fakulta EVOLUČNÍ ALGORITMY. Josef Tvrdík

UČEBNÍ TEXTY OSTRAVSKÉ UNIVERZITY Přírodovědecká fakulta

EVOLUČNÍ ALGORITMY Josef Tvrdík

OSTRAVSKÁ UNIVERZITA 2004

Obsah 1 Úvod

3

2 Problém globální optimalizace 2.1 Formulace problému globální optimalizace . . . . . . . . . . 2.2 Stochastické algoritmy pro globální optimalizaci . . . . . . .

4 4 5

3 Úvod do Matlabu 3.1 Co je Matlab? . . . . . . . . . . . . . . 3.2 Základní příkazy . . . . . . . . . . . . 3.3 Některé příkazové zkratky . . . . . . . 3.4 Znakové řetězce . . . . . . . . . . . . . 3.5 Nakreslení grafu . . . . . . . . . . . . . 3.6 Logické výrazy . . . . . . . . . . . . . 3.7 Řídicí příkazy, zápis algoritmu . . . . . 3.8 Doporučení pro vhodný výběr příkazů 3.9 3D grafika . . . . . . . . . . . . . . . . 3.10 Práce se soubory . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

8 8 8 10 13 13 14 15 18 20 22

4 Náhodné prohledávání

25

5 Simplexová metoda

28

6 Řízené náhodné prohledávání

31

7 Genetické algoritmy

35

8 Evoluční strategie

39

9 Diferenciální evoluce

43

10 Algoritmus SOMA

48

11 Evoluční prohledávání

53

12 Algoritmy se soutěžícími heuristikami

56

13 Experimentální ověřování a porovnávání algoritmů 13.1 Testovací funkce . . . . . . . . . . . . . . . . . . . . . . . . . 13.2 Experimentální ověřování algoritmů . . . . . . . . . . . . . .

62 62 67

14 Literatura

71

2

1

Úvod

Tento text je určen jako studijní opora povinně volitelného předmětu Evoluční algoritmy v kombinovaném (a případně i distančním) magisterském studiu na Přírodovědecké fakultě Ostravské university ve studijním oboru informační systémy. Cílem textu je poskytnout základy evolučních algoritmů, zejména těch, které jsou vhodné pro hledání globálního minima v souvislé oblasti a také seznámit studenty s principy programování v Matlabu, numerickým testováním stochastických algoritmů pro globální optimalizaci a vyhodnocením numerických experimentů. Každá kapitola začíná pokyny pro její studium. Tato část je vždy označena jako Průvodce studiem s ikonou na okraji stránky. V závěru každé kapitoly je rekapitulace nejdůležitějších pojmů z akuální kapitoly. Rekapitulace označena textem Shrnutí a ikonou na okraji. Oddíl Kontrolní otázky označený ikonou by vám měl pomoci zjistit, zda jste prostudovanou kapitolu pochopili a snad vyprovokuje i vaše další otázky, na které budete hledat odpověď. U některých kapitol je připometuto, že k dané kapitole budete řešit Korespondeční úlohu. Pokud nebudou korespondenční úlohy budou zadávány k vašemu kurzu jinak, je úspěšné vyřešení korespondenčních úloh z textu podmínkou pro získání zápočtu. Pojmy a důležité souvislosti k zapamatování jsou vyznačeny na okraji stránky textu ikonou.

3

2

Problém globální optimalizace

Průvodce studiem V této kapitole je zformulován problém globální optimalizace a stručně jsou zmíněny základní myšlenky stochastických algoritmů, které se používají pro nalezení globálního minima účelových funkcí. Počítejte asi se dvěma hodinami studia.

2.1

Formulace problému globální optimalizace

Úlohu nalezení globálního minima můžeme formulovat takto: Mějme účelovou funkci f : D → R,

D ⊆ Rd .

Pak x∗ = arg minx∈D f (x) je globální minimum. Globální minimum je tedy jeden bod nebo více bodů z D s nejmenší funkční hodnotou, formálně x∗ = arg min f (x) = {x∗ ∈ D : f (x∗ ) ≤ f (x) pro ∀x ∈ D} x∈D

Chceme-li nalézt globální maximum, pak jej nalezneme jako globální minimum funkce g(x) = −f (x). Problém nalezení globálního minima můžeme ilustrovat jednoduchým obrázkem:

4

Ze základního kurzu matematické analýzy známe postup, jak nalézt extrémy funkcí, u kterých existuje první a druhá derivace. Zdálo by se, že úloha nalezení globálního minima je velmi jednoduchá. Bohužel tomu tak není. Nalézt obecné řešení takto jednoduše formulovaného problému je obtížné, zvláště když účelová funkce je multimodální, není diferencovatelná, případně má další nepříjemné vlastnosti. Analýza problému globální optimalizace ukazuje, že neexistuje deterministický algoritmus řešící obecnou úlohu globální optimalizace (tj. nalezení dostatečně přesné aproximace x∗ ) v polynomiálním čase, tzn. problém globální optimalizace je NP-obtížný. Přitom globální optimalizace je úloha, kterou je nutno řešit v mnoha praktických problémech, mnohdy s velmi významným ekonomickým efektem, takže je nutné hledat algoritmy, které jsou pro řešení konkrétních problémů použitelné. Algoritmy pro řešení problému globální optimalizace se podrobně zabývá celá řada monografií, např. Torn a Žilinskas [28], Míka [20], Spall [26], kde je možné najít mnoho užitečných poznatků přesahujících rámec tohoto předmětu. V tomto kurzu se soustředíme na problémy, kdy hledáme globální minimum v souvislé oblasti Q D = ha1 , b1 i × ha2 , b2 i × . . . × had , bd i = di=1 hai , bi i, (2.1) ai < bi , i = 1, 2, . . . , d, (tzv. box constraints, protože oblast D je vymezena jako d-rozměrný kvádr) a účelovou funkci f (x) umíme vyhodnotit s požadovanou přesností v každém bodě x ∈ D. Pro úlohy řešené numericky na počítači nepředstavuje podmínka (2.1) žádné podstatné omezení, neboť hodnoty ai , bi jsou tak jako tak omezeny datovými typy užitými pro x a f (x).

2.2

Stochastické algoritmy pro globální optimalizaci

Nemožnost nalézt deterministický algoritmus obecně řešící úlohu globální optimalizace vedla k využití algoritmů stochastických, které sice nemohou garantovat nalezení řešení v konečném počtu kroků, ale často pomohou nalézt v přijatelném čase řešení prakticky použitelné. Stochastické algoritmy pro globální optimalizaci heuristicky prohledávají prostor D. Heuristikou rozumíme postup, ve kterém se využívá náhoda, intuice, analogie a zkušenost. Rozdíl mezi heuristikou a deterministickým algoritmem je v tom, že na rozdíl od deterministického algoritmu heuristika nezajišťuje nalezení řešení. Heuristiky jsou v praktickém životě zcela samozřejmě užívané postupy, jako příklady můžeme uvést hledání hub, lov ryb na udici, výběr partnera, pokus o výhru ve sportovním utkání nebo o složení zkoušky ve škole. 5

Většina stochastických algoritmů pro hledání globálního minima v sobě obsahuje zjevně či skrytě proces učení. Inspirace k užití heuristik jsou často odvozeny ze znalostí přírodních nebo sociálních procesů. Např. simulované žíhání je modelem pomalého ochlazování tuhého tělesa, tabu-search modeluje hledání předmětu tak, že v krátkodobé paměti si zapamatovává zakázané kroky vedoucí k již dříve projitým stavům. Popis těchto algoritmů najdete v knize Kvasničky, Pospíchala a Tiňa [15], případně v článcích [11, 12, 7, 8]. Podobné postupy učení najdeme snad ve všech známých stochastických algoritmech s výjimkou slepého náhodného prohledávání. V posledních desetiletích se s poměrným úspěchem pro hledání globálního minima funkcí užívají stochastické algoritmy zejména evolučního typu. Podrobný popis této problematiky naleznete v knihách Goldgergra [9], Michalewicze [18] nebo Bäcka [2], případně ve zprávě Bäcka a Schwefela [3], která je dostupná na Internetu. Rozvoj evolučních algoritmů je záležitostí posledních desetiletí a je podmíněn rozvojem počítačů a pokroky v informatice. Přelomovými pracemi zavádějící biologickou terminologii do modelů hledání globálního extrému jsou články Schwefela a Rechenberga ze šedesátých let minulého století o evoluční strategii, práce L. Fogela o evolučním programování - viz [2, 15] - a Hollandova kniha o genetických algoritmech [10]. Tím byl odstartován bouřlivý a rozsáhlý rozvoj evolučních algoritmů a jejich aplikací. Evoluční algoritmy jsou ve své podstatě jednoduchými modely Darwinovy evoluční teorie vývoje populací. Charakteristické pro ně je to, že pracují s populací a využívají tzv. evoluční operátory, zejména tyto: • selekce – nejsilnější jedinci z populace mají větší pravděpodobnost přežití a předání svých vlastností, • křížení (rekombinace) – dva nebo více jedinců z populace si vymění informace a vzniknou tak noví jedinci kombinující vlastnosti rodičů, • mutace – informace zakódovaná v jedinci může být náhodně modifikována. Evoluční algoritmy jsou heuristiky, které nějakým způsobem modifikují populaci tak, aby se její vlastnosti zlepšovaly. O některých třídách evolučních algoritmů je dokázáno, že nejlepší jedinci populace se skutečně přibližují ke globálnímu optimu. Evoluční algoritmy byly a jsou předmětem intenzivního výzkumu, počet a rozsah publikací z této oblasti je ohromující. Jedním z hlavních motivů jsou především aplikace v praktických problémech, které jinými metodami nejsou řešitelné. Rozsáhlé přehledy aplikací evolučních algoritmů jsou uvedeny např. v pracích [4] a [16]. Dalšími motivy jsou výzkum umělé inteligence a teorie učení. Hledají se nové inspirace např. v chovatelství a 6

šlechtitelství [22], modeluje se sexuální reprodukce v populacích [24] nebo chování mravenců [6]. Pokračuje výzkum sofistikovaných evolučních operátorů a přístupů [5, 19, 36, 37]. Tento rozvoj nezastavil ani tzv. ”No Free Lunch Theorem” [33, 34], který ukazuje, že nelze najít univerzálně nejlepší stochastický algoritmus pro globální optimalizaci.

Shrnutí

• Formulace problému globální optimalizace • Stochastický algoritmus • Heuristika Kontrolní otázky

1. Jak naleznete minimum diferencovatelné funkce analyticky? 2. Kdy jste ve svém životě užili heuristiku? Zkuste ji popsat.

7

3

Úvod do Matlabu

Průvodce studiem. Cílem této kapitoly je seznámit se s Matlabem tak, abyste v něm mohli pohodlně programovat evoluční algoritmy. Počítejte asi se čtyřmi hodinami studia, zejména praktických cvičení u počítače.

3.1

Co je Matlab?

Matlab je výpočetní prostředí pro maticové operace, technické výpočty a vizualizaci dat. Matlab už po řadu let vyvíjí firma MathWorks a má mnoho uživatelů v běžných operačních systémech (Windows, Unix) po celém světě. Je to jak interaktivní výpočetní systém, který interpretuje příkazy zadané v příkazovém okně, tak programovací jazyk vysoké úrovně, který umožňuje rychlou a jednoduchou implementaci algoritmů. Základním datovým typem jazyka je pole, což může být skalár, vektor, matice nebo dokonce více než dvourozměrné pole. Prvkem pole může být číselná hodnota v pohyblivé řádové čárce, dvojitá přesnost, rozsah od −10308 do 10308 (realmax=1.7977e+308, realmin=2.2251e-308), případně alfanumerický znak nebo řetězec znaků stejné délky. Logické hodnoty se vyjadřují číselně, false jako 0, true jako 1. Matlab kromě základních operací s takovým polem a spousty vestavěných dalších funkcí nabízí i programovací příkazy pro vytváření podprogramů (funkcí) a příkazy pro řízení průběhu výpočtu (podmíněný příkaz, switch, cykly) známé z běžných programovacích jazyků. Výhodou Matlabu je i velmi kvalitní dokumentace, která je přístupná online prostřednictvím Helpu a dovoluje jak rychlé vyhledávání témat, zejména funkcí Matlabu, tak rychlou a názornou instruktáž k základním operacím. Pro seznámení s ovládáním Matlabu a jeho základními možnostmi je velmi užitečné spustit Help a projít Getting Started.

3.2

Základní příkazy

Pro užívání Matlabu je důležité zvládnout základní operace s vektory a maticemi a také naučit se „myslet vektorověÿ, abychom byli připraveni efektivně využít možnosti, které máme v Matlabu k dispozici. Nejdříve si ukážeme pár jednoduchých příkladů. Řádkový vektor a vytvoříme příkazem, ve kterém prvky vektoru oddělené mezerami nebo čárkou zadáme v hranatých závorkách. Po zadání příkazu se tento vektor také vypíše v příkazovém okně: >> a=[3 2 5] a = 3

2

5 8

Pokud bychom chtěli výpis potlačit, příkaz ukončíme středníkem. Tím jsme vytvořili v pracovní oblasti (Workspace) objekt s názvem odpovídajícím užitému identifikátoru, t.j. a, který můžeme dále užívat v dalších příkazech. Pozor! Matlab v identifikátorech rozlišuje malá a velká písmena, takže a a A nebo TEMP_1, temp_1 a Temp_1 jsou různé objekty. Existující objekt můžete užít např. jako vstupní parametr funkce: >> size(a) ans = 1

3

Funkce size vrací rozměry matice, první údaj je počet řádků, druhý počet sloupců. Příkazem >> b=a+1 b = 4

3

6

jsme přičetli zadanou skalární hodnotu (t.j. 1) ke každému prvku vektoru. Matici vytvoříme např. příkazem, ve kterém jsou řádkové vektory odděleny středníkem >> B=[a; b; b-5.5] B = 3.0000 2.0000 4.0000 3.0000 -1.5000 -2.5000

5.0000 6.0000 0.5000

Funkce length vrací maximální hodnotu z vektoru spočítaného funkcí size: >> length(B) ans = 3 Prvky matice nebo vektoru jsou přístupné přes jejich indexy, takže např. prvek na druhém řádku a v prvním sloupci můžeme přepsat: >> B(2,1)=0 B = 3.0000 0 -1.5000

2.0000 3.0000 -2.5000

5.0000 6.0000 0.5000

9

S vektory a maticemi můžeme snadno provádět operace známé z lineární algebry (transpozice, inverze matice, pokud je možná, atd.) Jejich zápis je velmi podobný zápisu obvyklém v lineární algebře: >> a=a’ a = 3 2 5 >> det(B) % determinant ans = 54 >> C=B^-1 % inverse matice B, stejný výsledek užitím inv(B) C = 0.3056 -0.2500 -0.0556 -0.1667 0.1667 -0.3333 0.0833 0.0833 0.1667 >> I=B*C I~= 1.0000 0 0.0000

3.3

0.0000 1.0000 -0.0000

0 0 1.0000

Některé příkazové zkratky

Pro vytvoření vektorů a matic jsou užitečné další příkazy, které zkracují práci. Pokud tvoří prvky vektoru sekvenci, pak ji můžeme zadat zkráceně zadáním dolní a horní hranice (pak je implicitně krok roven jedné) nebo požadovaný krok explicitně zadat: >> a=1:5 a = 1 >> b=10:-2:1 b = 10

2

3

4

5

8

6

4

2

10

Matice můžeme vytvářet pomocí příkazů matice (n × p), samé nuly matice (n × p), samé jedničky jednotková matice (n × n) matice (n × p), náhodná čísla z rovnoměrného rozdělení na [0, 1) matice (n × p), náhodná čísla z normovaného normálního rozdělení

zeros(n,p) ones(n,p) eye(n) rand(n,p) randn(n,p)

Zadáme-li jen jeden parametr, pak vytvořená matice je čtvercová. Příklady: >> zeros(1,3) ans = 0 >> eye(4) ans = 1 0 0 0 >> rand(3) ans = 0.4447 0.6154 0.7919

0

0 1 0 0

0

0 0 1 0

0 0 0 1

0.9218 0.7382 0.1763

0.4057 0.9355 0.9169

>> 10+2*randn(3) ans = 7.1181 11.3800 11.1423 11.6312 9.2002 11.4238

12.5805 11.3372 12.3817

Pro práci s vektory a maticemi je důležité prázdné pole [], které má rozměr (0×0). Zařazením prázdného prvku (matice) vypouštíme odpovídající prvky příslušného pole: a = 1 >> a(2)=[] a = 1

2

3

4

3

4

5

5

11

Řádek nebo sloupec matice můžeme přepsat vektorem >> X=ones(3,4) X = 1 1 1 1 1 1 >> X(1,:)=a X = 1 3 1 1 1 1

1 1 1

1 1 1

4 1 1

5 1 1

nebo prázdným polem, tj. vypustit >> X(:,3)=[] X = 1 3 1 1 1 1

5 1 1

Pole můžeme i zvětšovat (přidávat nové prvky). Pokud by v novém poli nebyla nějaká hodnota dosud přiřazena, automaticky se dodají nuly, jak uvidíme v následujícím příkladu: a = 1 >> a(5)=-1 a = 1

2

3

2

3

0

-1

Zvláštní význam má proměnná end obsahující aktuální nejvyšší hodnotu indexu. Pak např. přepsání posledního prvku vektoru a hodnotou předposledního prvku zapíšeme takto: >> a(end)=a(end-1) a = 1 2

3

0

0

Kromě dosud zmíněných aritmetických operátorů jsou užitečné operace na odpovídajících prvcích dvou polí stejných rozměrů, tj. operátory začínající tečkou, .*, ./, .^. Např. umocnit každý prvek vektoru a odpovídajícím prvkem téhož vektoru zapíšeme takto a dostaneme výsledek: >> a.^a ans =

1

4

27

1

1 12

Pokud výsledek aritmetické operace „přetečeÿ rozsah datového typu, je výsledkem hodnota Inf, např. >> realmax*1.2 ans = Inf >> -realmax*1.2 ans = -Inf Pokud výsledek aritmetické operace nelze vyhodnotit (výrazy Inf*0 nebo Inf/Inf ap.), má hodnotu NaN (not a number), např. >> -realmax*1.2/(realmax*1.2) ans = NaN

3.4

Znakové řetězce

Znakový řetězec v Matlabu je řádkový vektor, jehož prvky jsou jednotlivé znaky, např. >> r2=’abcdA1B2’ r2 =abcdA1B2 length(r2) ans = 8 Pokud chceme přímo porovnávat dva řetězce, musí mít shodnou délku. Shodné řetězce musí mít stejnou délku a shodovat se na všech pozicích. Pokud potřebujete porovnávat řetězce nestejných délek, užijte vhodnou funkci, např. strcmp, viz Help, kde najdete i návody pro další operace se znakovými řetězci.

3.5

Nakreslení grafu

V Matlabu je velmi snadné kreslení grafů funkcí. Zde si ukážeme jen jednoduchý příklad 2D grafu: >> x=0:pi/400:6*pi; >> y=exp(-0.2*x).*sin(x); >> plot(x,y) Zadáním těchto tří příkazů dostaneme v okně Figure(1) následující graf:

13

Podrobnější popis příkazu plot, jeho dalších nepovinných parametrů, možností interaktivní práce s grafem, export obrázku do různých formátů atd. nalezneme v helpu. Některé postupy vytváření 3D grafů ukážeme na příkladech později.

3.6

Logické výrazy

Logické konstanty false a true se vyjařují numericky jako 0 a 1, při vyhodnocování výrazu je jakákoliv číselná hodnota různá od nuly považována za true. Logické spojky not, and, or se zapisují jako ∼, &&, ||. Zdvojením && a || je zajištěno zrychlené vyhodnocování logického výrazu (zleva). Pro exkluzivní nebo je k dispozici funkce xor. Relační operátory se zapisují očekávaným způsobem, tj. ==, ∼=, <, <=, >, >=, rovnítko musí být zdvojeno, aby bylo odlišitelné od přiřazovacího příkazu. V Matlabu máme k dispozici i kvantifikátory any a all. Je-li jejich parametrem vektor, vracejí skalár, tzn. logickou hodnotu, zda aspoň jeden, resp. všechny prvky vektoru mají logickou hodnotu true. Např. a = 1 2 >> all(a) ans = 1

3

14

Pokud je jejich parametrem matice, vracejí řádkový vektor pro jednotlivé sloupce matice, např. pro A~= 0 1 0 0 >> any(A) ans = 0 >> all(A) ans = 0

1 1 1

1

0

1

Poznámka: Podobně fungují i některé numerické funkce, např. mean, std, median, min, max, sum, prod, takže můžeme velmi snadno spočítat řádkový vektor sloupcových průměrů, směrodatných odchylek atd. Pomocí operátorů any a all můžeme snad zjistit, zda aspoň jeden či všechny prvky matice nebo vektoru splňují nějakou logickou podmínku, např. zda všechny prvky matice jsou nezáporné: B = -1 3 3 -1 0 3 >> all(all(B>=0)) ans = 0 nebo >> all(B(:)<=0) ans = 0 Nalezení indexů prvků vektoru splňujících logickou podmínku umožňuje funkce find. Např. záporné prvky vektoru b vypustíme snadno takto b = 0.2621 -0.0435 >> b(find(b<0))=[] b = 0.2621 0.3214

3.7

-0.4815

0.3214

-0.0553

Řídicí příkazy, zápis algoritmu

V Matlabu máme k dispozici běžné příkazy pro zápis algoritmu, známé z jiných programovacích jazyků. Sekvence příkazů je buď na následujícím řádku nebo lze příkazy oddělit středníkem. Je-li příkaz ukončen středníkem, vyhodnocený výraz se nevypisuje na displej. Pokud má zápis příkazu pokračovat na dalším řádku, je předcházen třemi tečkami. Identifikátor muže 15

mít až 31 alfanumerických znaků nebo podtržítek, začíná písmenem, velká a malá písmena se rozlišují. Jako identifikátory neužívejte Matlabem užívané názvy konstant a proměnných jako např. pi, eps, ans, Inf, NaN, realmax, realmin, neboť tím byste si přepsali jejich správné nastavené hodnoty. Vše za znakem procento je komentář. Podmíněný příkazy má tvar if podmínka příkaz1 elseif podmínka2 příkaz2 else příkaz3 end Části elseif a else jsou nepovinné, podmíněný příkaz může být zapsán i na jednom řádku. Přepínač má tvar switch(výraz) case value1 příkaz1 case value2 příkaz2 .. . case valuen příkazn otherwise nějaký příkaz end

Část otherwise je nepovinná. Cykly lze zapsat buď jako while podmínka příkaz end nebo for identifikátor = vektor příkaz end Všude, kde je možno zapsat příkaz, může to být i sekvence více příkazů. Např. v cyklu typu for proběhne desetkrát přiřazení vektoru náhodných čísel do i−tého řádku matice: 16

for i=1:10 P(i,:)=2+3*rand(1,5); end Sekvence příkazů je možné zapsat do textových souborů, jejichž jméno má příponu „.mÿ. Celou sekvenci příkazů pak vyvoláme jménem souboru bez přípony. Pokud tato sekvence nemá žádné vstupní parametry, říkáme takovému souboru skript (dávka). Při jejím vyvolaní je pak sekvence příkazů postupně interpretována. Velmi užitečným prostředkem je zapsat sekvenci příkazů jako funkci. Tak máme k dispozici postup pro modulární strukturu programů v Matlabu. Syntaxi zápisu funkce můžeme ve stručnosti vyjádřit schématem function [seznam výstupních parametrů] =jmeno_funkce(seznam vstupních parametrů) Pokud funkce vrací jen jednu hodnotu (jednu proměnnou, může to být ale i vektor nebo matice), výstupní parametr není (nemusí být) v hranatých závorkách. Proměnné zavedené uvnitř funkce jsou lokální, takže nepřepisují hodnoty stejně nazvaných proměnných ve volajícím modulu. Zvolené jméno funkce musí být shodné se jménem souboru, ve kterém je zdrojový text funkce uložen (bez přípony „.mÿ). Užívání funkcí je i výhodné s ohledem na délku výpočtu, neboť při prvním volání proběhne kompilace funkce a další volání je rychlejší než opakovaná interpretace jednotlivých příkazů. Jako příklad funkce si ukážeme náhodný výběr k hodnot (bez opakování) z množiny {1, 2, . . . , N }, k ≤ N . Pokud by bylo N = 49 a k = 7, je to procedura losování Sportky. Příklad pečlivě prostudujte, neboť ilustruje několik užitečných operací s vektory v Matlabu. % random sample, k of N without repetition % function result=nahvyb(N,k); opora=1:N; nahv=[]; for i=1:k index=1+fix(rand(1)*length(opora)); nahv(end+1)=opora(index); opora(index)=[]; end result=nahv; Uvedený způsob implementace náhodného výběru je jen jedním z mnoha možných. Je to vlastně velmi přímočarý model fyzického losování z osudí. Jiný způsob implementace může využít funkci randperm. Vyzkoušejte to a posuďte, který postup je efektivnější. 17

Při volání funkce nemusí být vždy zadán plný počet vstupních i výstupních parametrů. Pro ošetření takového volání jsou určeny funkce nargin a nargout (bez parametrů), které vracejí aktuální počet příslušných parametrů. Vstupním parametrem může být i jméno funkce (řetězec), takže příkazem y=feval(fnname,x); vyhodnotíme v bodě x tu funkci, jejíž jméno je aktuálně uloženo jako řetězec v proměnné fnname. Samozřejmě funkce tohoto jména musí být dostupná, např. jako m-soubor.

3.8

Doporučení pro vhodný výběr příkazů

Při zápisu algoritmů v Matlabu se vyplatí „uvažovat vektorověÿ, neboť řadu operací je možno v Matlabu zapsat efektivněji než kopírováním postupů známých z jazyků jako je Pascal nebo C. Především je nutno velmi uvážlivě užívat cyklů, zejména vložených. Různou časovou náročnost ukazují následující příklady výpočtu součtu druhých mocnin prvků vektoru a. Časově nejnáročnější je „klasický postupÿ načítání druhých mocnin jednotlivých prvků v cyklu: N=10000 a=rand(50,1) tic for i=1:N s=0; for j=1:length(a) s=s+a(j)^2; end end toc elapsed_time = 1.2700 Využitím aritmetické operace pro všechny prvky vektoru (zde .* nebo .^2) dosáhneme téměř osminásobného zkrácení potřebného času: tic for i=1:N s=sum(a.*a); end toc elapsed_time =

0.1600 18

A využijeme-li skalární součin, potřebný čas ještě o trochu zkrátíme: tic for i=1:N s=a’*a; end toc elapsed_time =

0.1100

19

3.9

3D grafika

Kreslení grafů funkcí v 3D si ukážeme jen na jednoduchých příkladech. V prvním příkladu chceme nakreslit graf funkce f (x, y) = exp(−(x2 + y 2 )) pro hodnoty x ∈ [−1, 2; 1.2] a y ∈ [−2; 2]. Nejdříve vytvoříme vektory hodnot x-ové i y-ové souřadnice (pro síť vhodné hustoty): x=-1.2:0.05:1.2 y=-2:0.1:2 Pak příkazem meshgrid vytvoříme matice X a Y typu (length(y) × length(x)), v tomto případě tedy typu (41 × 49), ve kterých odpovídající prvky představují všechny dvojice hodnot prvků vektorů x a y. Pak už jen vyhodnotíme z-tovou souřadnici (matice Z typu (41 × 49)) [X Y]=meshgrid(x,y) Z=exp(-X.^2-Y.^2) mesh(X,Y,Z) % stejného výsledku dosáhneme i zadaním příkazu % ve tvaru mesh(x,y,Z) a příkazem mesh nakreslíme obrázek (tzv. drátový graf), který se objeví v okně Figure 1.

V tomto okně máme pak dostupné nástroje pro interaktivní práci s obrázkem, např. popis a editaci os grafu, vkládání textů a grafických symbolů, rotace obrázku, uložení do souboru, export do jiného formátu atd. Podrobnosti 20

ponecháváme na čtenáři, který může potřebné informace získat v HELPu, kde také se dozví i o dalších možnostech 3D grafiky v Matlabu. O trochu složitější je nakreslení trojrozměrného grafu funkce tehdy, když matici z-tových souřadnic nemůžeme spočítat přímo, ale máme jen funkci vyhodnocující z-tovou souřadnici v bodě zadaném souřadnicemi x a y. Příkladem je Ackleyho funkce, která je často užívána v testování algoritmů pro globální optimalizaci a kterou v Matlabu můžeme zapsat takto: function z=ackley(x) % Ackley’s function, <-30,30>, glob. minimum 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); z=a+b+20+exp(1); Síť souřadnic v rovině xy vytvoříme podobně jako v předchozím příkladu x=-3:0.05:3; y=-3:0.05:3; [X Y]=meshgrid(x,y); ale výpočet z-tových souřadnic umíme pouze v bodech zadaných jako dvouprvkový vektor. Pak nezbývá než spočítat matici Z po jednotlivých prvcích for i=1:length(y) for j=1:length(x) Z(i,j)=ackley([X(i,j) Y(i,j)]); end end Výpočet je možno trochu zrychlit, když z matic X a Y vytvoříme vektory, výsledný vektor zz předem inicializujeme na požadovanou délku a pak z něho příkazem reshape vytvoříme matici požadovanou pro kreslení 3D grafu: xx=X(:); yy=Y(:); zz=zeros(length(xx),1); for i=1:length(xx) zz(i)=ackley([xx(i) yy(i)]); end Z=reshape(zz,length(y),length(x)); Pak teprve můžeme nakreslit graf následujícími příkazy: 21

meshc(X,Y,Z) axis([-3 3 -3 3 -1 4]) xlabel(’x’) ylabel(’y’) zlabel(’z’)

3.10

Práce se soubory

V Matlabu máme řadu možností, jak importovat data ze souborů a exportovat data do souborů různých formátů. Podrobnosti nalezneme v HELPu pod heslem file formats a odkazech na další položky. Zde uvedeme jen některé příklady práce se soubory, které se budou užitečné při interaktivní práci s Matlabem, při načítání vstupních dat z textových souborů a při ukládání výsledků do textové tabulky se sloupci s konstantní šířkou. Při interaktivní práci s Matlabem je někdy užitečné si uložit aktuální stav pracovního prostoru (workspace), tj. všechny dosud inicializované proměnné, abych v příštím sezení mohli pokračovat od současného stavu. To je možné volbou save Workspace as z menu File. Soubor se pak uloží pod zadaným jménem s příponou .MAT. Tento soubor pak otevřeme z menu File ve volbě Open. Matici (6 × 2) z textového souboru, jehož jméno je v řetězci (proměnná) fdatname a soubor má následující obsah 22

2.138E0 3.421E0 3.597E0 4.340E0 4.882E0 5.660E0

1.309E0 1.471E0 1.490E0 1.565E0 1.611E0 1.680E0

načteme snadno do proměnné (matice) X příkazem dlmread X=dlmread(fdatname,’ ’); Druhý parametr znamená specifikaci oddělovače (delimiter) datových položek, v tomto případě je to mezera (alespoň jedna mezera znamená konec číselné položky, u jiných oddělovačů jako je tabulátor nebo středník to je jinak, tam každý výskyt oddělovače znamená novou položku). Podobně snadno lze načítat data i z jiných formátů, např. tabulky z Excelu ap. Zapisování dat do textového souboru s pevnou šířkou sloupců zařídíme sekvencí příkazů: Otevření souboru příkazem fid=fopen(fname,’a’); kde první parametr je jméno souboru (proměnná nebo konstanta typu řetězec), druhý parametr specifikuje, jakým způsobem má být soubor otevřen, v tomto případě hodnota parametru ’a’ znamená append, tj. přidání dalších záznamů, pokud soubor už existuje. Pak následují příkazy pro zápis jednotlivých položek řádku souboru. Vždy seznam položek, které mají být zapsány, je přecházen specifikací formátu, poslední příkaz ukončuje řádek ve výstupním souboru. fprintf(fid,’ %-10s %4.0f %11.3e %11.3e’,mod_name,d,a(1),b(1)); fprintf(fid,’ %7.2f’ , alfa ); fprintf(fid,’ %8.0f’, evals); fprintf(fid,’ %15.7e’,RSS); for i=1:d fprintf(fid,’ %14.5e’, beta(i)); end fprintf(fid,’%1s\n’,’ ’); Další záznam do souboru pak zapíšeme voláním stejné sekvence příkazů po změně obsahu výstupních proměnných, obvykle tedy tato sekvence je umístěna uvnitř nějakého cyklu. Na závěr soubor zavřeme příkazem fclose(fid). Uvedené příklady jen ilustrují základní možnosti práce se soubory, o dalších operacích se soubory včetně práce s grafickými formáty se potřebné informace naleznou v HELPu. 23

Shrnutí

• Datové typy v Matlabu, prázdné pole • Inicializace vektorů a matic • Skripty a funkce v Matlabu • 3D grafy Kontrolní otázky

1. Jak budete generovat náhodné číslo z rovnoměrného rozdělení z intervalu [a, b)? 2. Kdy může být jméno funkce vstupním parametrem jiné funkce? 3. Jak generovat náhodný výběr k hodnot z [1, 2, . . . , N ] pomocí funkce randperm?

24

4

Náhodné prohledávání

Průvodce studiem V této kapitole se seznámíte s nejjednodušším stochastickým algoritmem pro hledání globálniho minima. Kromě toho uvidíte i snadnost zápisu takových algoritmů v Matlabu a seznámíte se s pojmem konvergence stochastického algoritmu. Počítejte asi se dvěma hodinami studia. Náhodné prohledávání se také někdy nazývá slepý algoritmus, neboť v něm se Qdopakovaně generuje náhodné řešení (nový bod y v souvislé oblasti D = i=1 hai , bi i z rovnoměrného spojitého rozdělení) a zapamatovává se tehdy, když je lepší než předtím nalezemé řešení. Náhodné prohledávání je příklad velmi jednoduchého stochastického algoritmu pro hledání globálního minima využívající velmi prostou heuristiku, kterou bychom mohli slovně vyjádřit jako „zkus bez rozmyslu cokolivÿ. Slepý algoritmus lze zapsat v Matlabu na pá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 the minimal function value found % x minimum found by 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 Kromě jednoduchosti je výhodou tohoto algoritmu také to, že lze dokázat jeho konvergenci. Platí, že s rostoucím počtem iterací se nalezené řešení x přibližuje globálnímu minimu x∗ : lim P (k (x − x∗ ) k<  | t) = 1  > 0,

t→∞

(4.1)

kde k (x − x∗ ) k je vzdálenost nalezeného řešení x od globálního minima x∗ (norma vektoru x − x∗ ) a P (k (x − x∗ ) k<  | t) je pravděpodobnost jevu (k (x − x∗ ) k< ) za podmínky, že bylo provedeno t iterací. 25

K tvrzení vyjářenému rovnicí (4.1) dojdeme následující úvahou: Všechny body splňující podmínku k (x − x∗ ) k<  jsou vnitřními body d-rozměrné koule D o poloměru ). Jelikož  > 0, je „objemÿ této koule kladný, přesněji míra této množiny λ(D ) > 0. Označme dále míru množiny D jako λ(D). Generujeme-li bod y ∈ D z rovnoměrného rozdělení, pak pravděpodobnost jevu A ≡ {y ∈ D } λ(D ) P (A) = λ(D) a pravděpodobnost jevu opačného je ¯ =1− P (A)

λ(D ) λ(D)

Pravděpodobnost, že jev A¯ nastane v t po sobě jdoucích opakováních je t  λ(D )  P (A¯ | t) = 1 − . λ(D)   ) Jelikož 0 < 1 − λ(D < 1, je potom λ(D) lim P (A¯ | t) =

t→∞



λ(D ) 1− λ(D)

t =0

a tedy rovnice(4.1) platí. Bohužel rovnice (4.1) vypovídá pouze limitně a neposkytuje nám žádnou informaci o kvalitě nalezeného řešení po konečném počtu iterací. Jak uvidíme později, s podobným problémem se budeme setkávat i u dalších stochastických algoritmů. Teoretická konvergence podle rovnice (4.1) je z praktického hlediska velmi slabé tvrzení a pro posouzení konvergence stochastických algoritmů a porovnání jejich efektivity jsou většinou důležitější experimentální výsledky na testovacích funkcích než teoretický důkaz konvergence ve formě rovnice (4.1).

26

Shrnutí

• Generování bodu z rovnoměrného rozdělení na D • Konvergence stochastického algoritmu Kontrolní otázky

1. Jakou heuristiku užívá slepý algoritmus pro generování nového bodu? 2. Je ve slepém náhodném prohledávání modelován nějaký proces učení? Korespondenční úloha Zjistěte experimentálně, jaká je pro první DeJongovu funkci (viz kapitola 13) průměrná hodnota t potřebná k tomu, aby f (x) < 0.05

27

5

Simplexová metoda

Průvodce studiem V této kapitole se seznámíte s algoritmem zvaným simplexová metoda. Je to příklad algoritmu, který hledá minimum funkce ve spojité oblasti a nevyužívá derivací funkce. Řada dalších stochastických algoritmů využívá některé postupy, které byly zavedeny v simplexové metodě. Proto věnujte této metodě pozornost, zejména jasnému porozumění pojmu simplex a reflexe simplexu. Počítejte zhruba se dvěma až třemi hodinami studia. Simplexová metoda je velmi jednoduchý a populární algoritmus pro hledání globálního minima. Původní verzi algoritmu publikovali Nelder a Mead v roce 1965 [23]. Od té doby bylo navrženo několik dalších modifikací simplexové metody. Zde si ukážeme její základní myšlenky. Pro účelovou funkci D ⊂ Rd .

f : D → R,

Q hledáme její globální minimimum v souvislé oblasti D = di=1 hai , bi i, a účelovou funkci f (x) umíme vyhodnotit s požadovanou přesností v každém bodě x ∈ D. Funkce nemusí být diferencovatelná. Klíčovým pojmem algoritmu je tzv. simplex S, což je množina nekomplanárních d + 1 bodů v D, S = {x1 , x2 , . . . , xd+1 } V simplexu nalezneme body s nejvyšší a nejnižší funkční hodnotou xH = arg max f (x) x∈S

xL = arg min f (x) x∈S

a dále spočítáme těžiště, tj. průměr d bodů, které zůstanou v simplexu po odstranění bodu xH , tedy 1 X g= ( x − xH ) d x∈S Pomocí tohoto simplexu hledáme nový bod y, který v případě, že f (y) < f (xH ) nahradí v simplexu dosud nejhorší bod xH a pokračuje se znovu. Pro hledání nového bodu y se užívá reflexe. Reflexe znamená překlopení bodu xH přes těžiště g podle vztahu y = g + (g − xH ) = 2g − xH . Graficky je reflexe znázorněna na následujícím obrázku.

28

(5.1)

Pokud nový bod y získaný reflexí nesplňuje podmínku f (y) < f (xH ), užije se redukce, což je vlastně takové „smrštěníÿ simplexu směrem k bodu xL , že vzdálenosti bodů simplexu budou poloviční. Formálně to zapíšeme 1 x ← (x + xL ), x ∈ S, x 6= xL . 2 Na následujícím obrázku vidíme grafické znázornění redukce. Body původního simlexu označené kroužky se posunou k bodu xL do pozic označených křížky.

Nejjednodušší variantu algoritmu simplexové metody můžeme zapsat následujícím způsobem: 29

generuj simplex S, tj. d + 1 bodů náhodně v D repeat y := reflexe(S); (y ∈ D) if f (y) < f (xH ) then xH := y else redukce(S); until podmínka ukončení;

Verbálně můžeme simplexovou metodu jako jakýsi „píďalkovitýÿ pohyb simplexu v prohledávaném prostoru D směrem k oblasti s nižšími funkčními hodnotami. Mezi variantami simplexové metody existují i takové, kde simplex na počátku není generován náhodně, ale počáteční body simplexu se záměrně volí tak, aby zcela jistě splňovali podmínku nekomplanarity (afinní nezávislosti), což volně řečeno znamená, že body simplexu neleží v podprostoru s menší dimenzí, než má prohledávaný prostor D, např. při d = 2 tři body tvořící simplex neleží na přímce. Sofistikovaná verze simplexové metody je užita ve funkci fminserch, která je standardní součástí Matlabu. Prostudujte help k této funkci a její zdrojový text (soubor fminserch.m), ve kterém najdete nejen popis implementované verze simplexové metody, ale také užitečnou inspiraci pro styl programování v Matlabu.

Shrnutí

• Simplex • Reflexe v simplexu • Redukce simplexu Kontrolní otázky

1. Naprogramujte v Matlabu jednoduchou verzi simplexové metody popsanou v této kapitole a porovnejte počet vyhodnocení účelové funkce se slepým náhodným prohledáváním na úloze z předcházející kapitoly. Korespondenční úloha Napište seznam nejdúležitějších odlišností algoritmu užitém ve funkci fminserch od jednoduché simplexové metody popsané v této kapitole. 30

6

Řízené náhodné prohledávání

Průvodce studiem V této kapitole se seznámíte s algoritmem, který pracuje s populací. Algoritmus řízeného náhodného prohledávání je příkladem jednoduchého a přitom efektivního algoritmu pro hledání globálního minima. Pochopení tohoto algoritmu je dobrým odrazovým můstkem ke složitějím evolučním algoritmům. Počítejte asi se dvěma hodinami studia a několika hodinami věnovanými implementaci algoritmu v Matlabu a jeho ověření na testovacích funkcích. Řízené náhodné prohledávání (controlled random search, CRS) je příkladem velmi jednoduchého a přitom efektivního stochastického algoritmu pro hledání minima v souvislé oblasti. Algorimus CRS navrhl Price [25] v sedmdesátých létech minulého století. Pracuje se s populací N bodů v prohledávaném prostoru D a z nich se heuristicky generuje nový bod y, který může být zařazen do populace místo dosud nejhoršího bodu. Počet vygenerovaných bodů populace N je větší než dimenze d prohledávaného prostoru D. Jako heuristiku pro generování nového bodu y užíval Price reflexi simplexu, která je známa ze simplexové metody. Z populace se vybere náhodně d + 1 nekomplanárních bodů tvořících simplex. Pokud v novém bodu y je funkční hodnota f (y) menší než je v nejhorším bodu populace, pak je tento nejhorší bod nahrazen bodem y. Reflexi v simplexu můžeme vyjádřit vztahem y = 2g − x,

(6.1)

kde x je bod simplexu s největší hodnotou účelové funkce a g je těžiště zbývajících d bodů simplexu. Nahrazením nejhoršího bodu populace novým bodem y dosahujeme toho, že populace se koncentruje v okolí dosud nalezeného bodu s nejmenší funkční hodnotou. Algoritmus mužeme jednoduše zapsat následujícím způsobem: generuj populaci P, tj. N bodů náhodně v D repeat najdi xworst ∈ P takové, že f (xworst ) ≥ f (x), repeat y := reflexe simplexu, y ∈ D until f (y) < f (xworst ); xworst := y; until podmínka ukončení;

x∈P

Reflexe simplexu podle rovnice (6.1) samozřejmě není jediná možnost, jak v algoritmu CRS generovat nový bod. Price později navrhl postup, kdy nový bod y se generuje podle (6.1), ale do simplexu je vždy zařazen nejlepší bod 31

populace xmin s funkční hodnotou fmin , fmin ≤ f (x), x ∈ P a zbývajících d bodů simplexu se pak vybere náhodně z ostatních bodů populace. Další možnost je znáhodněná reflexe, která byla navržena v 90. létech. Znáhodněná reflexe v simplexu je popsána vztahem y = g + U (g − xH ),

(6.2)

U je náhodná veličina vhodného rozdělení. Graficky je taková reflexe znázorněna na následujícím obrázku.

Pokud bychom užili reflexi podle (6.1), byla by vzdálenost bodů g, y byla stejná jako vzdálenost bodů x, g, tj. U v rovnici (6.2) je rovno jedné. Při užití znáhodněné reflexe je tato vzdálenost určována náhodně podle zvoleného rozdělení náhodné veličiny U . Na testovacích úlohách bylo ověřováno několik rozdělení náhodné veličiny U . Osvědčilo se rovnoměrné rozdělení na [0, α), kde α je vstupní parametr algoritmu. Podle empirických výsledků testů nejrychleji algoritmus konvergoval na většině testovacích funkcí při hodnotách α ≈ 4. Střední hodnota je E U = α/2. Za pozornost stojí užití rovnoměrného rozdělení náhodné veličiny U na [s, α − s), kde α > 0 a s ∈ (0, α/2) jsou vstupní parametry heuristiky. Střední hodnota je opět E U = α/2. Pro další výzkum připadají v úvahu i jiná rozdělení produkující nezáporné hodnoty náhodné veličiny U , např. lognormální rozdělení, případně i jiné heuristiky než reflexe. 32

Reflexe podle rovnic (6.1) nebo (6.2) nezaručuje, že nově vygenerovaný bod y bude v prohledávaném prostoru D. Pokud nastane situace, že y 6∈ D, jsou dvě možnosti, jak postupovat dále. Buď generujeme Qd nový bod opakovaně tak dlouho, než se podaří vygenerovat y ∈ D = i=1 hai , bi i nebo použít výpočetně úspornější postup tzv. zrcadlení (perturbace), kdy ty souřadnice yi 6∈ hai , bi i překlopíme dovnitř prohledávaného prostoru D kolem příslušné strany d-rozměrného kvádru D. Toto zrcadlení znázorňuje následující obrázek, jeho algoritmus je zapsán jako funkce zrcad v Matlabu.

% zrcadleni, perturbation y into function result=zrcad(y,a,b) zrc=find(yb); for i=zrc while (y(i)b(i)) if y(i)>b(i) y(i)=2*b(i)-y(i); elseif y(i)
33

Při stochastickém prohledávání podmínku ukončení formuluje uživatel. Obvykle je možné užít podmínku ukončení formulovanou tak, že nejlepší a nejhorší bod populace se ve funkčních hodnotách liší jen málo, tedy f (xworst ) − f (xmin ) ≤ ε,

(6.3)

kde ε > 0 je vstupní parametr algoritmu, případně podmínka ukončení může být formulována jako disjunkce uvedené relace (6.3) a dosažení maximálního dovoleného počtu vyhodnocení účelové funkce.

Shrnutí

• Prohledávaný prostor D, způsob zadání D • Populace bodů náhodně vygenerovaných v D • Simplex náhodně vybraný z populace • Reflexe a znáhodněná reflexe • Podmínka ukončení Kontrolní otázky

1. Jak určíte souřadnice těžiště n bodů v D? 2. Co je to zrcadlení a kdy se užívá? 3. Jak zformulovat podmínku ukončení? 4. Naprogramujte v Matlabu algoritmus CRS

34

7

Genetické algoritmy

Průvodce studiem Tato stručná kapitola má posloužit k základní orientaci v principech genetických algoritmů a upozornit na postupy, které jsou inspirativní pro jiné evoluční algoritmy. Počítejte asi se dvěma hodinami studia. Genetické algoritmy od svého vzniku v 70. létech [10] představují mohutný zdroj inspirace pro další evoluční algoritmy, pro zkoumání umělé inteligence a pro rozvoj soft computing. O genetických algoritmech existuje řada specializovaných monografií, viz např. Goldberger [9], Michalewicz [18] nebo Bäck [2]. Velmi dobře a srozumitelně napsaný úvod do genetických algoritmů najdete v dostupné knize Kvasničky, Pospíchala a Tiňa [15]. Zde se omezíme jen na stručné vysvětlení základních principů těchto algoritmů. Základní myšlenkou genetických algoritmů je analogie s evolučními procesy probíhajícími v biologických systémech. Podle Darwinovy teorie přirozeného výběru přežívají jen nejlépe přizpůsobení jedinci populace. Mírou přizpůsobení je tzv. „fitnessÿ jedince. V biologii je fitness chápána jako relativní schopnost přežití a reprodukce genotypu jedince. Biologická evoluce je změna obsahu genetické informace populace v průběhu mnoha generací směrem k vyšším hodnotám fitness. Jedinci s vyšší fitness mají větší pravděpodobnost přežití a větší pravděpodobnost reprodukce svých genů do generace potomků. Kromě reprodukce se v populačním vývoji uplatňuje i tzv. mutace, což je náhodná změna genetické informace některých jedinců v populaci. V genetických algoritmech je fitness kladné číslo přiřazené genetické informaci jedince. Tato genetická informace jedince (chromozom) se obvykle vyjadřuje bitovým řetězcem. Populaci jedinců pak modelujeme jako populaci chromozomů a každému chromozomu (bitovému řetězci) umíme přiřadit hodnotu účelové funkce a podle vhodného předpisu umíme přiřadit i jeho fitness, podrobněji viz [15]. Obvykle chromozomem je binární vektor konstantní délky k α = (α1 , α2 , . . . , αk ) ∈ {0, 1}k a populace velikosti N je potom množina takových chromozomů P = {α1 , α2 , . . . , αN } Pokud hledáme globální minimum účelové funkce f , pak fitness F je nějaké zobrazení, které vyhovuje následující podmínce f (αi ) ≤ f (αj ) ⇒ F (αi ) ≥ F (αj ) > 0, i, j = 1, 2, . . . , N Nejjednodušší způsob, jak vyhovět této podmínce, je lineární zobrazení funkčních hodnot na fitness, kdy hodnotu fitness určíme podle vztahu F (α) =

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

kde fmin a fmax je nejmenší a největší 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é číslo, např.  = 0.01. Pak vztah pro fitness má tvar F (α) =

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

Pak je vhodné hodnoty fitness renormalizovat, aby jejich součet byl roven jedné. Renormalizovaná fitness i-tého jedince je F (αi ) F 0 (αi ) = PN j=1 F (αj ) a tato hodnota pak udává i pravděpodobnost účasti jedince na reprodukci jeho genetické informace, tj. vytvoření potomka. Výběr (selekci) jedince s pravděpodobností F 0 (αi ) ukazuje následující algoritmus, vstupem je vektor fitness všech jedinců populace, hodnoty fitness nemusí být renormalizovány. 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(1)=fits(1); for i=2:h cp(i)=cp(i-1)+fits(i); end cp=cp/ss; res=1+fix(sum(cp l1 a v bitových řetězcích se vymění úseky αl1 , . . . , αl2 ↔ βl1 , . . . , βl2 36

Mutace je většinou v genetických algoritmech implementována jako změna hodnoty náhodně vybraného bitu v chromozomu, tj. hodnota 0 se změní na 1, hodnota 1 na 0. Mutací je zajištěno, aby mohla v populaci vzniknout genetická informace, která v něm v předcházejících generací nebyla, či obnovení genetické informace ztracené v průběhu dosavadního vývoje. V důsledku toho pak algoritmus může uniknout z oblasti lokálního minima, ke kterému by směřoval, pokud bychom užívali pouze křížení. Naopak, pokud bychom užívali pouze mutaci, genetický algoritmus by se choval podobně jako slepé náhodné prohledávání. Lze ukázat, že genetický algoritmus konverguje ke globálnímu extrému, viz [15]. S genetickými algoritmy a tímto důkazem jste se také setkali v předmětu Neuronové sítě, takže zde se jím nebudeme zabývat. Genetické algoritmy na rozdíl od ostatní algoritmů uváděných v tomto textu reprezentují jedince jako bitový řetězec, nikoliv jako vektor hodnot v pohyblivé řádové čárce. To je často výhodou při řešení tzv. diskrétních úloh globální optimalizace (kombinatorické úlohy jako je problém obchodního cestujícího), ale přínáší to komplikace v problémech, prohledávaná oblast D možných řešení je spojitá. Pokud bychom úseky chromozomu interpretovali jako celá čísla (datový typ integer), pak se potýkáme s obtíží, že sousední číselné hodnoty se jsou reprezentovány řetězci, které se mohou podstaně lišit, dokonce i ve všech bitech, např. dekadická hodnota 7 je binárně (0111), následující číselná hodnota 8 je (1000). Tuto nepříjemnost je nutno v genetických algoritmech nějak ošetřit, pokud mají být úspěšně používány i pro řešení spojitých problémů globální optimalizace. Je to možno řešit buď použitím Grayova kódu [15], kdy řetězce reprezentující sousední celočíselné hodnoty se liší jen v jednom bitu, ale přináší to další časový nárok na konverzi do Grayova kódu a zpět nebo se závádějí varianty genetických algoritmů s reprezentací jedince jako vektoru číselných hodnot v pohyblivé čárce. Pak je ovšem nutné užívat jiné operace křížení a mutace než ty, které jsou uvedeny výše.

37

Shrnutí

• Evoluce v biologických systémech • Genetická informace, chromozom, fitness • Operátory selekce, křížení a mutace v genetických algoritmech Kontrolní otázky

1. Proč jako fitness nelze užívat přímo hodnotu účelové funkce? 2. Jaká je role mutace v genetických algoritmech?

38

8

Evoluční strategie

Průvodce studiem V této kapitole jsou vysvětleny základy evoluční strategie, která patří k nejstarším a stále velmi populárním a inspirujícím evolučním algoritmům. Počítejte asi se třemi hodinami studia. Původní nejjednodušší verzi algoritmu navrhli v šedesátých létech Schwefel a Rechenberg. Základní myšlenka je podobná slepému náhodnému prohledávání, ale rozdíl je v tom, že nový bod y se generuje jako mutace bodu x tak, že jednotlivé složky vektoru x se změní přičtením hodnot normálně rozdělených náhodných veličin y = x + u,

u ∼ N (0, σ 2 I),

(8.1)

tj. u = (U1 , U2 , . . . , Ud ) je náhodný vektor, jehož každý prvek je náhodná veličina Ui ∼ N (0, σ 2 ), i = 1, 2, . . . , d a tyto veličiny jsou navzájem nezávislé. Algoritmus lze zapsat v Matlabu takto: function [x,fx] = es1p1const(fn_name,a,b,t,sigma) % input parameters: % fn_name function to be minized (M file) % a, b row vectors, limits of search space % t iteration number % 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); fy= feval(fn_name,y); if fy < fx x=y; fx=fy; end end Evoluční terminologií může být algoritmus popsán tak, že z rodičovské generace velikosti 1 vzniká generace potomků rovněž velikosti 1, potomek vzniká mutací z jednoho rodiče podle rov. (8.1) a operárorem selekce je výběr lepšího jedince z dvojice (tzv. turnajový výběr). Zajímavou otázkou je to, jak volit hodnoty směrodatných odchylek pro mutaci. V rov. (8.1) jsou všechny σj , j = 1, 2, . . . , d shodné a konstantní. To 39

samozřejmě není nutné, každá dimenze může mít svou hodnotu směrodatné odchylky, tedy budeme pracovat s vektorem směrodatných odchylek σ = (σ1 , σ2 , . . . , σd ) a navíc, vektor směrodatných odchylek nemusí být v každém iteračním kroku stejný, ale hodnoty jeho prvků se mohou adaptovat podle průběhu procesu vyhledávání. Rechenberg studoval vliv velikosti mutace při generování potomka, tj. hodnot σ na rychlost konvergence algoritmu a na dvou jednoduchých funkcích odvodil tzv. pravidlo jedné pětiny úspěšnosti. Empiricky se pak byla ověřena užitečnost tohoto pravidla i pro jiné funkce. Hodnoty směrodatných odchylek se v i-té iteraci hledání upravují podle pravidla vyjádřeného rovnicí  1  c1 σi−1 když ϕ(n) < 5 σi = (8.2) c σ když ϕ(n) > 15  2 i−1 σi−1 jinak c1 < 1 a c2 > 1 jsou vstupní parametry, kterými se řídí zmenšování či zvětšování hodnot směrodatných podle relativní četnosti úspěchu ϕ(n) v předcházejících n krocích. Úspěchem se rozumí to, že f (y) < f (x). Počet kroků n je vstupní parametr. Obvykle se volí hodnoty c1 = 0.82 a c2 = 1/c1 = 1.22. Pod vlivem rozvoje jiných evolučních algoritmů bylo dalším krokem v evoluční strategii zavedení populací velikostí větších než jedna. Označíme-li velikost rodičovské populace M a velikost populace potomků L, L > M , pak můžeme uvažovat o dvou variantách algoritmu evoluční strategie, v literatuře obvykle označovaných jako ES(M+L) a ES(M,L): • ES(M+L) - generace potomků je vytvořena M jedinci s nejmenšími funkčními hodnotami ze všech M +L jedinců jak rodičovské populace, tak populace potomků. V této variantě se dodržuje tzv. elitismus, tj. nejlepší jedinec vždy přežívá a přechází do nové generace. • ES(M,L) - generace potomků je vytvořena M jedinci s nejmenšími funkčními hodnotami z L jedinců populace potomků. Rodičovské generace tedy je kompletně nahrazena nejlepšími jedinci z populace potomků. V literatuře (viz např. [15]) se doporučuje, aby velikost populace potomků L byla volena několikanásobně větší než velikost rodičovské populace M . Také se uvádí, že varianta ES(M,L) obvykle konverguje pomaleji než ES(M+L), ale s menší tendencí ukončit prohledávání v lokálním minimu. Algoritmus ES(M+L) v Matlabu je uveden na následující stránce.

40

function [x,fx,func_evals] =... es_mpl(fn_name,a,b,M,k,lfront,my_eps,max_evals) d=length(a); sigma=(b-a)/3; % initial values of sigma fr=rand(1,lfront) < 0.2; % initial values in last lfront steps L=M*k; % k>1, integer P=ones(M,d+1); Q=ones(L,d+1); 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)); 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 0.2 sigma=sigma*1.22; end end end PQ=[P; Q]; PQ=sortrows(PQ,d+1); P=PQ(1:M,:); P(1,d+1) end % main loop - end x=P(1,1:d); fx=P(1,d+1);

41

Autoři algoritmu evoluční strategie v minulosti museli čelit výtkám, že evoluční strategie z osvědčených evolučních operátorů evoluční strategie vůbec nevyužívá křížení. Proto byly navrženy pokročilejší varianty evoluční strategie, ve kterých je křížení obsaženo. Základní idea takového křížení spočívá v tom, že u každého jedince v populaci je kromě souřadnic v prostoru D i jeho vektor směrodatných odchylek a vektor směrodatných odchylek potomka se pak generuje křížením s náhodně vybraným jiným jedincem. Máme-li dva rodiče a, b s vektory směrodatných odchylek σ a , σ b , pak vektor směrodatných odchylek jejich potomka je určován např. jako σ=

σa + σb , 2

což je tzv. křížení průměrem nebo podle nějaké obdobného jen o trochu komplikovanějšího pravidla pro křížení. Pro funkce, u kterých je jejich globální minimum v protáhlém údolí, jehož směr není rovnoběžný se žádnou dimenzí prostoru D, se užívá ještě sofistikovanější varianta evoluční strategie, v níž u každého jedince v populaci se kromě souřadnic v prostoru D a jeho vektoru směrodatných odchylek uchovávají i mimodiagonální prvky kovarianční matice (tzv. směrové úhly, viz např. [15]) a křížení probíhá nejen na vektorech směrodatných odchylek, ale na celé kovarianční matici. Takové varianty evoluční strategie jsou však nejenom implemetačně podstatně náročnější, ale také mají větší paměťové nároky i větší časovou spotřebu na jeden iterační krok. Shrnutí

• mutace v evoluční strategii • adaptace hodnot směrodatných odchylek • ES(M,L) a ES(M+L) Kontrolní otázky

1. Co je pravidlo jedné pětiny a k čemu se užívá? 2. Porovnejte výhody a nevýhody variant ES(M,L) a ES(M+L). 3. Jak je pravidlo jedné pětiny implementováno ve variantě ES(M+L) v Matlabu a kterým parametrem řídíme délku sledovaného počtu předchozích iterací?

42

9

Diferenciální evoluce

Průvodce studiem V této kapitole se seznámíte s algoritmem diferenciální evoluce. Tento algoritmus byl navržen nedávno a poprvé publikován v roce 1995. Během několika let se stal velmi populární a je často aplikován na problémy hledání globálního minima. Je to příklad jednoduchého heuristického hledání užívajícího evoluční evoluční operátory. Počítejte asi se třemi hodinami studia. Diferenciální evoluce (DE) je postup k heuristickému hledání minima multimodálních funkcí, který navrhli R. Storn a K. Price [27] koncem 90. let. Algoritmus DE dosáhl značné popularity, je často aplikován a dále rozvíjen. Bibliografický přehled prací o diferenciální evoluci naleznete na webových stránkách J. Lampinena. Experimentální výsledky i zkušenosti z četných aplikací ukazují, že často konverguje rychleji než jiné stochastické algoritmy pro globální optimalizaci. Algoritmus diferenciální evoluce vytváří novou populaci Q tak, že postupně pro každý bod xi ze staré populace P vytvoří jeho potenciálního konkurenta y a do nové populace z této dvojice zařadí bod s nižší funkční hodnotou. Algoritmus můžeme zapsat takto: generuj P (N bodů náhodně v D) repeat for i := 1 to N do generuj vektor u vytvoř vektor y křížením u a xi if f (y) < f (xi ) then do Q zařaď y else do Q zařaď xi endfor P := Q until podmínka ukončení Je několik způsobů, jak generovat nový bod u. Zde uvedeme dva nejčastěji užívané. Postup označovaný RAND generuje bod u ze tří bodů ze staré populace podle vztahu u = r1 + F (r2 − r3 ), (9.1) r1 , r2 , r3 jsou navzájem různé body náhodně vybrané z populace P různé od aktuálního bodu xi , F > 0 je vstupní parametr. Generování bodu u podle rovnice 9.1 je graficky znázorněno následujícím obrázkem.

43

Postup označovaný BEST využívá nejlepší bod z populace P a čtyři další náhodně vybrané body podle vztahu u = xbest + F (r1 + r2 − r3 − r4 ),

(9.2)

kde r1 , r2 , r3 , r4 jsou navzájem různé body náhodně vybrané z populace P různé od aktuálního bodu xi i od bodu xbest s nejnižší funkční hodnotou v P . F > 0 je opět vstupní parametr. Nový vektor y vznikne křížením vektoru u a vektoru xi tak, že kterýkoli jeho prvek xij je nahrazen hodnotou uj s pravděpodobností C. Pokud žádné xij nebylo přepsáno hodnotou uj nebo při volbě C = 0, nahrazuje se jeden náhodně vybraný prvek vektoru xi :  uj když Rj ≤ C nebo j = I yj = (9.3) xij když Rj > C a j 6= I , kde I je náhodně vybrané celé celé číslo z {1, 2, . . . , d}, Rj ∈ (0, 1) jsou voleny náhodně a nezávisle pro každé j a C ∈ [0, 1] je vstupní parametr. V diferenciální evoluci generování u podle rovnic (9.1) nebo (9.2) představuje mutaci, křížení probíhá podle (9.3), selekcí je výběr lepšího bodu z dvojice y a xi do nové populace Q. V diferenciální evoluci jsou tedy přítomny všechny tři nejdůležitější evoluční operace.

44

Ali a Törn nedávno navrhli určovat hodnotu parametru F v každém iteračním kroku podle adaptivního pravidla ( max(Fmin , 1 − | ffmax |) if | ffmax |<1 min min F = (9.4) fmin max(Fmin , 1 − | fmax |) jinak, kde fmin , fmax jsou minimální a maximální funkční hodnoty v populaci P a Fmin je vstupní parametr, který zabezpečuje, aby bylo F ∈ [Fmin , 1). Autoři doporučují volit hodnotu Fmin ≥ 0.45. Předpokládá se, že tento způsob výpočtu F udržuje prohledávání diverzifikované v počátečním stadiu a intenzivnější v pozdější fázi prohledávání, což má zvyšovat spolehlivost hledání i rychlost konvergence. Výhodou algoritmu diferenciální evoluce je jeho jednoduchost a výpočetně nenáročné generování nového bodu y, které lze navíc velmi efektivně implementovat. Nevýhodou je poměrně velká citlivost algoritmu na nastavení hodnot vstupních parametrů F a C. Storn a Price doporučují volit hodnoty N = 10d, F = 0.8 a C = 0.5 a pak tyto hodnoty modifikovat podle empirické zkušenosti z pozorovaného průběhu hledání. To je ovšem doporučení dosti vágní. K tomu, aby bylo užitečné, je nutná zkušenost a dobrá intuice řešitele problému. V testovacích úlohách v článku [27] sami užívají pro různé úlohy velmi odlišné hodnoty těchto parametrů, a to 0.5 ≤ F ≤ 1 a 0 ≤ C ≤ 1. Také velikost populace N volí často menší než doporučovanou hodnotu N = 10d. Porovnáme-li diferenciální evoluci s algoritmem CRS, vidíme, že v diferenciální evoluci se nenahrazuje nejhorší bod v populaci, ale pouze horší bod ve dvojici, což zabezpečuje větší diversitu bodů nové populace po delší období, takže při stejné velikosti populace většinou DE ve srovnání s CRS má menší tendenci ukončit prohledávání v lokálním minimu, ale za to platíme pomalejší konvergencí při stejné podmínce ukončení. Jednoduchost implementace DE s generováním nového bodu postupem RAND ilustruje následující zdrojový text v Matlabu.

45

function [x_star,fn_star,func_evals]=... der_smp(fn_name, a, b, N, my_eps, max_evals, F, CR) % 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))); 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=derand(P,F,CR,i); y=zrcad(y,a,b); % perturbation fy=feval(fn_name,y); 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;

46

function y=derand(P,F,C,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); v=r1+F*(r2-r3); change=find(rand(1,d)
Shrnutí

• Mutace jako přičtení rozdílu vektorů • Křížení výměnou náhodně vybraných prvků vektoru • Výběr lepšího z dvojice (selekce turnajem) Kontrolní otázky

1. Jaké jsou hlavní odlišnosti řízeného náhodného výběru a diferenciální evoluce? 2. Proč musí být F 6= 0? Co by se stalo, kdybychom zadali F < 0? 3. Co způsobí zadání C = 0? 47

10

Algoritmus SOMA

Průvodce studiem V této kapitole se stručně seznámíte s algoritmem SOMA, který modeluje chování skupiny mající společný cíl. Počítejte asi se třemi až čtyřmi hodinami studia. Algoritmus SOMA byl publikován nedávno [36] a dosti podrobně, ale nepříliš přehledně je popsán v české knize jednoho z autorů tohoto algoritmu [37]. SOMA je zkratka dlouhého jména algoritmu (Self-Organizing Migration Algorithm), které vystihuje jeho základní principy. Na algoritmus můžeme nahlížet jako na jednoduchý model lovící smečky. Ve smečce je N jedinců a ti se pohybují (migrují) po území D tak, že každý jedinec skáče směrem k vůdci smečky, prověří místo doskoku a zapamatuje si nejlepší místo do dalšího migračního kola. Vůdcem je většinou jedinec s nejlepší hodnotou účelové funkce v celé populaci (smečce) a rovněž kvalita míst (bodů v D) navštívených jedincem při jeho skocích za vůdcem se posuzuje hodnotou účelové funkce. Pohyb migrujícího jedince nemusí probíhat ve všech dimenzích prostoru D, ale jen v některých dimenzích, tedy v prostoru s dimenzí menší než d. Migraci jedince k vůdci ukazuje následující obrázek. Jedinec z výchozího bodu r0 skáče skoky velikosti δ směrem k vůdci, případně ho může i o trochu přeskočit. Pohyb jedince je buď ve směru naznačeném plnou šipkou, tedy v prostoru dimenze d nebo jen v prodprostoru menší dimenze, tedy v některém ze směrů vyznačených tečkovaně.

Kromě vstupních parametrů obvyklých u všech evolučních algoritmů, tj. 48

specifikace problému (funkce, prohledávaný prostor, podmínka ukončení) a velikost populace N má algoritmus SOMA ještě tři vstupní parametry popisující uvedený pohyb jedince za vůdcem. Jsou to tyto parametry (jejich označení ponecháme shodné s tím, které užívá ve své knize Zelinka [37]: • mass – relativní velikost přeskočení vůdce 0

k xl − r0 k , mass = k xl − r0 k 0

kde k xl − r0 k je norma vektoru (xl − r0 ), podobně k xl − r0 k 0 je norma vektoru (xl − r0 ). Tyto normy vektoru jsou vlastně délky 0 úseček s koncovými body r0 , xl , resp. r0 , xl . Geometrický význam je také zřejmý z obrázku. Doporučeny jsou hodnoty mass ∈ [1.1; 3]. • step – určuje velikost skoků δ, tj. velikost směrového vektoru skoku podle následujícího vztahu δ = step ∗ (xl − r0 ), tedy čím je hodnota parametru step menší, tím podrobněji je prostor na cestě k vůdci prozkoumáván. Doporučené hodnoty jsou step ∈ [0.11; mass], autoři algoritmu SOMA rovněž doporučují volit hodnotu parametru step tak, aby 1/step nebylo celé číslo, tzn. aby jedinec na své cestě za vůdcem nenavštívil zbytečně místo obsazené vůdcem. • prt – parametr pro určení směru pohybu jedince za vůdcem, prt ∈ [0; 1]. Je to pravděpodobnost výběru dimenze prostoru D, ve které bude pohyb jedince probíhat. Pro směr pohybu se vyberou náhodně ty dimenze, pro které Ui < prt, i = 1, 2, . . . , d. Hodnoty (U1 , U2 , . . . , Ud ) tvoří náhodný vektor velikosti d, jehož složky jsou nezávislé náhodné veličiny rovnoměrně rozdělené na intervalu [0, 1). Při imlementaci algoritmu se tyto hodnoty generují funkcí rand. Pokud by žádná složka Ui nesplňovala podmínku Ui < prt (např. při volbě prt = 0), vybere se jedna dimenze náhodně z množiny {1, 2, . . . , d}. Směr pohybu, tzn. ty dimenze, ve kterých probíhají skoky, se volí vždycky pro každý pohyb jedince k vůdci znovu, zatímco velikost skoků, tedy δ je stejná pro celý běh algoritmu SOMA.

49

Autoři rozlišují tři základní varianty algoritmu SOMA: • all-to-one – základní verze algoritmu SOMA, kdy v každém migračním kole všichni jedinci kromě vůdce putují za vůdcem. Tato verze je popsána v předchozím textu a její implementace v Matlabu je uvedena na konci této kapitoly. • all-to-all – v každém migračním kole se vůdcem stávají postupně všichni jedinci a ostatní putují za nimi. Populace se aktualizuje nově nalezenými body až po skončení migračního kola, tzn. po dokončení putování všech ke všem. • all-to-all-adaptive – podobně jako u předchozí varianty v každém migračním kole se vůdcem stávají postupně všichni jedinci a ostatní putují za nimi, ale jedinec v populaci je nahražen bezprostředně po ukončení jeho migrace za vůdcem (pokud při této migraci byl nalezen bod s lepší funkční hodnotou), takže ostatní jedinci v průběhu migračního kola už putují k tomuto novému vůdci.

50

% SOMA all to 1 function [x_star,fn_star,func_evals]=... soma_to1(fn_name,a,b,N,my_eps,max_evals,mass,step,prt) % input parameters: % obligatory: % 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 % output: % fn_star the minimal function value found by GCRS % x_star 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)); end % 0-th generation initialized P=sortrows(P,d+1); pocet_kroku=fix(mass/step); func_evals=N; while ((P(N,d+1)-P(1,d+1)) > my_eps) ... &(func_evals
51

% 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); d=length(r0); tochange=find(rand(d,1)>prt); if length(tochange)==0 % at least one is changed tochange=1+d*fix(rand(1)); 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); C(i,:)=[r, fr]; end [fmin, indmin]=min(C(:,d+1)); yfy=C(indmin,:); % row with min. fy is returned

Shrnutí

• Model pohybu jedince za vůdcem a parametry mass, step, prt • Varianta all-to-one • Varianta all-to-all • Varianta all-to-all-adaptive Kontrolní otázky

1. Proč parametr step má být volen tak, aby 1/step nebylo celé číslo? 2. Jak se liší varianta all-to-all od all-to-all-adaptive Korespondenční úloha Naprogramujte v Matlabu variantu all-to-all (jde o poměrně snadnou modifikaci varianty all-to-one uvedené v učebním textu), porovnejte alespoň na dvou testovacích funkcích (viz kap. 13) s variantou all-to-one. Pošlete stručnou zprávu o výsledcích testů a m soubor obsahující ověřenou variantu all-to-all. 52

11

Evoluční prohledávání

Průvodce studiem V této kapitole se budeme zabývat algoritmem vycházejícím ze řízeného náhodného prohledávání a obohaceným o další evoluční operátory. Počítejte asi se dvěma hodinami studia. Evoluční prohledávání (Evolutionary Search) je modifikací algoritmu řízeného náhodného prohledávání (CRS), v níž je posílena role evolučních oprerátorů. Nevýhodou řízeného náhodného prohledávání je to, že nejhorší jedinec je z populace odstraněn v každém iteračním kroku. Tím algoritmus CRS ztrácí schopnost uniknout z oblasti lokálního minima a řízené náhodné prohledávání má někdy tendenci k tzv. přečasné konvergenci, tj. ukončení v lokálním minimu (prodrobněji viz kapitola 13.2). V evolučním prohledávání se pracuje se dvěma populacemi P a Q, obě populace mají stejnou velikost N . Ze staré populace P se vybere M nejlepších bodů, které přejdou do nové populace. Tím se provede selekce nejlepších jedinců pro přežití. Za reprodukci můžeme považovat generování „potomkůÿ nějakou heuristikou, např. reflexí v simplexu, který se náhodně vybírá ze staré populace P . Potomek je do nové populace Q zařazen tehdy, je-li lepší než nejhorší jedinec ve staré populaci P . Po doplnění nové Q na velikost N se pro další generaci stane starou populací, v té se najde nejvyšší funkční hodnota, která bude sloužit jako kritérium pro zařazeni nového jedince do populace Q. Pak pokud je splněna podmínka pro mutaci, tak se jednomu z vybraných bodů staré populace P nádhodně změní hodnoty jeho prvků, tedy souřadnice v prostoru D. Tím se pak tento zmutovaný jedinec účastní v této generaci reprodukce (např. může být vybrán do simplexu) a podílet se tak na případném úniku z oblasti lokálního minima. Způsob selekce užitý v tomto algoritmu je podobný jako v evoluční strategii typu ES(N + (N − M )).

53

Algoritmus evolučního prohledávání můžeme zapsat takto: generate P (an old population of N points taken at random from D) find the worst point in P, xworst , with the highest value of f repeat copy M best points of P into new population Q, 1 ≤ M < N repeat repeat generate a new trial point y by some heuristics until f (y) < f (xworst ) insert the next trial point into Q until Q is completed to N points replace P by Q find the worst point in P, xworst , with the highest value of f if condition for mutation is true then mutate a point in P until stopping condition is true Vstupní parametry algoritmu evolučního prohledávání jsou stejné jako u CRS a navíc ještě počet jedinců M přecházejících do nové populace a pravděpodobnost mutace pmut . Připomeňme, že algoritmus CRS je speciálním případem algoritmu evolučního prohledávání. Zvolíme-li hodnotu vstupních parametrů M = N − 1 a pravděpodobnost mutace rovnou nule (podmínka pro mutaci má pak vždy hodnotu ”false”), pak evoluční prohledávání je ekvivalentní algoritmu CRS. Algoritmus evolučního prohledávání byl navržen roku 1996 [29] a osvědčil se jako velmi spolehlivý při odhadu parametrů nelineárních regresních modelů v úlohách, kde klasické deterministické algoritmy selhávájí a také algoritmus CRS nebyl stoprocentně spolehlivý [14]. Cenou za spolehlivost je však v případě evolučního prohledávání pomalejší konvergence, tj. k dosažení podmínky ukončení je potřeba větší počet vyhodnocení účelové funkce, většinou zhruba dvakrát až třikrát větši v porovnání s CRS. Pro evoluční prohledávání lze najít podmínky, za kterých algoritmus konverguje s pravděpodobností rovnou jedné. Z výsledků [21] vyplývá, že evoluční prohledávání je konvergentní, když (a) nejlepší bod staré populace vždy přechází do nové populace (tzn. M ≥ 1, tzv. elitismus). (b) mutace zaručuje kladnou pravděpodobnost vygenerování nového bodu y ∈ S, kde S je libovolná otevřená podmnožina D mající Lebesgueovu míru λ(S) > 0. To je splněno, když jako mutace je užito náhodné generování bodu y ze spojitého rozdělení, jehož sdružená hustota fy (y) > 0 pro ∀y ∈ D. Této podmínce vyhovuje např. rovnoměrné rozdělení užívané ve slepém náhodném prohledávání. 54

(c) pravděpodobnosti P mutace pmut,k v každé z k generací hledání tvoří divergentní řadu, ( ∞ k=1 pmut,k = ∞). Tuto podmínku lze snadno splnit tím, že pravděpodobnost mutace volíme konstantní a kladnou, pmut > 0.

Shrnutí

• Evoluční prohledávání jako zobecnění algoritmu CRS • Evoluční operátory a jejich užití v algoritmu evolučního prohledávání • Podmínky konvergence algoritmu evolučního prohledávání Kontrolní otázky

1. Jak je v evolučním prohledávání realizována selekce a jak křížení? 2. Jak implementovat mutaci, aby byl algoritmus konvergentní? 3. Zdůvodněte podrobněji, proč je v porovnání s CRS evoluční prohledávání spolehlivější a pomalejší. Korespondenční úloha Naprogramujte v Matlabu evoluční prohledávání a porovnejte je alespoň na dvou testovacích funkcích (viz kap. 13) s algoritmem CRS. Pošlete stručnou zprávu o výsledcích testů a m soubory obsahující ověřované algoritmy evolučního prohledávání a CRS.

55

12

Algoritmy se soutěžícími heuristikami

Průvodce studiem Tato kapitola je věnována modifikacím algoritmů evolučního a řízeného náhodného prohledávání, ve kterých se heuristiky užívané pro generování nového bodu střídají. Dosahuje se tak adaptibility algoritmů na řešený problém a tím i rozšíření oblasti úspěšných aplikací. Počítejte asi s pěti hodinami studia. V této kapitole se seznámíme s dalším zobecněním algoritmů evolučního a řízeného náhodného prohledávání, ve kterých se heuristiky užívané pro generování nového bodu střídají a pravděpodobnost výběru konkrétní heuristiky je úměrná její dosavadní úspěšnosti. Principy této soutěže vysvětlíme na zobecněném algoritmu řízeného náhodného prohledávání. Algoritmus 1. Zobecněný CRS generuj populaci P, tj. N bodů náhodně v D; repeat najdi xmax ∈ P, f (xmax ) ≥ f (x), x ∈ P; repeat užij nějakou heuristiku k vygenerování nového bodu y ∈ D; until f (y) < f (xmax ); xmax := y; until podmínka ukončení; Nový bod y se v klasické variantě algoritmu CRS generuje reflexí v simplexu podle (6.1). Samozřejmě je však možné ke generování bodu y užít i jinou heuristiku, jak už bylo zmíněno v kapitole o algoritmu CRS, např. znáhodněnou reflexi podle (6.2). Bod y lze však generovat jakoukoli jinou heuristikou. Odtud už je jen krůček k užítí více heuristik v rámci jednoho algorimu a tyto heuristiky střídat. A pokud v tomto střídání heuristik budou preferovány ty, které byly v dosavadním průběhu hledání úspěšnější, je naděje, že algoritmus bude konvergovat rychleji. V takovém algoritmu heuristiky soutěží o to, aby byly vybrány pro generování nového bodu. Mějme tedy k dispozici h heuristik a v každém kroku ke generování nového bodu vybíráme náhodně i-tou heuristiku s pravděpodobností qi , i = 1, 2, . . . , h. Pravděpodobnosti qi měníme v závislosti na úspěšnosti i-té heuristiky v průběhu vyhledávacího procesu. Heuristiku považujeme za úspěšnou, když generuje nový bod y takový, že f (y) < f (xmax )). Různým hodnocením této úspěšnosti dostaneme různé varianty algoritmu se soutěžícími heuristikami.

56

Pokud ni je dosavadní počet úspěchů i-té heuristiky, pravděpodobnost qi je úměrná tomuto počtu úspěchů qi = Ph

ni + n0

j=1 (nj

+ n0 )

,

(12.1)

kde n0 > 0 je vstupní parametr algoritmu. Nastavením n0 ≥ 1 zabezpečíme, že jeden náhodný úspěch heuristiky nevyvolá příliš velkou změnu v hodnotě qi . Algoritmus užívající k hodnocení úspěšnosti heuristik (12.1) je jednou z možných variant. Jinou možností, jak ohodnotit úspěšnost heuristiky, je vážit úspěšnost relativní změnou v hodnotě funkce. Váha wi se určí jako wi =

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

(12.2)

Hodnoty wi jsou v intervalu (0, 1) a pravděpodobnost qi se pak vyhodnotí jako Wi + w0 , (12.3) qi = Ph j=1 (Wj + w0 ) kde Wi je součet vah wi v předcházejícím hledání a w0 > 0 je vstupní parametr algoritmu. Aby se zabránilo potlačení možnosti výběru některé z heuristik, lze zadat vstupní parametr δ a klesne-li kterákoli hodnota qi pod hodnotu δ, jsou pravděpodobnosti výběru heuristik nastaveny na jejich počáteční hodnoty qi = 1/h. Další prostor pro tvorbu různých variant algoritmu je volba heuristik, které spolu budou soutěžit. Příkladem heuristik jsou reflexe a znáhodněná reflexe v simplexu popsané v kapitole o algortimu CRS. Jiný příkladem jsou heuristiky vycházející z diferenciální evoluce, kdy bod u se generuje jako v diferecíální evoluci a nový vektor y vznikne křížením vektoru u a náhodně vybraného vektoru x. Heuristiky inspirované evoluční strategií generují nový bod y podle následujícího pravidla yi = xi + Y, i = 1, 2, . . . , d, kde xi je i-tá souřadnice náhodně vybraného bodu x z populace a Y je náhodná veličina, Y ∼ N (0, σi2 ). Zatímco v evoluční strategii je hodnota parametru σi2 adaptována v průběhu procesu vyhledávání podle tzv. pravidla 1/5 úspěchu, zde může být aktuální hodnota tohoto parametru odvozována z aktuální populace. Např. v heuristice es-pop je hodnota směrodatné odchylky σi určována z celé populace jako σi = k(max xi − min xi ) + ε, x∈P

x∈P

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

operátory max, min znamenají největší, resp. nejmenší hodnotu i-té souřadnice v aktuální populaci P , ε > 0 (v implementaci můžete užít např. 57

ε = 1 × 10−4 ) zabezpečuje, že hodnota σi je kladná, k je vstupní řídící parametr heuristiky, k > 0. V heuristice es-2pts je směrodatná odchylka σi určována ze dvou náhodně vybraných bod˚ u r, s σi = k |ri − si | + ε,

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

k je opět vstupní parametr heuristiky, ostatní symboly mají stejný význam jako výše. Heuristikou může být samozřejmě také náhodné generování bodu y ze spojitého rovnoměrného rozdělení na D, které se užívá ve slepém náhodném prohledávání. První varianta evolučního algoritmu se soutěžícími heuristikami byla publikována v práci [30], v práci [31] je podrobnější analýza tohto algoritmu. Výsledky ukázaly, že pro sestavení vhodné sady heuristik pro soutěž je užitečné znát vlastnosti jednotlivých heuristik a jejich role v procesu hledání. Zajímavé je, že evoluční algoritmus se soutěžícími heuristikami může být konvergentní, i když v něm není zahrnuta explicitně mutace (pmut = 0). To znamená, že konvergentní může být i Algoritmus 1, uvedený na začátku této kapitoly. Z výsledků [21] vyplývá, že takový evoluční algoritmus se soutěžícími heuristikami je konvergentní, když (a) mezi soutěžícími heuristikami je alespoň jedna heuristika zaručující kladnou pravděpodobnost vygenerování nového bodu y ∈ S, kde S je libovolná otevřená podmnožina D mající Lebesgueovu míru λ(S) > 0. (b) pravděpodobnosti výběru takové P heuristiky qk v každém z k kroků hledání tvoří divergentní řadu ( ∞ k=1 qk = ∞). (c) nejlepší bod staré populace vždy přechází do nové populace (tzv. elitismus). Podmínku (a) lze zabezpečit zařazením heuristiky, která náhodně generuje bod y ze spojitého rozdělení, jehož sdružená hustota fy (y) > 0 pro ∀y ∈ D. Této podmínce vyhovuje např. rovnoměrné rozdělení užívané ve slepém náhodném prohledávání. Podmínky (b) a (c) jsou zajištěny volbou vstupních parametrů algoritmu, tj. u evolučního prohledávání M ≥ 1 (u algoritmu CRS je podmínka (b) splněna automaticky) a kladné limitní hodnoty δ pro změnu dosažených hodnot pravděpodobnosti výběru soutěžících heuristik na počáteční hodnoty v průběhu vyhledávacího procesu. Pro usnadnění vaší implementace algoritmů se soutěžícími heuristikami v Matlabu následuje text dvou funkcí, které můžete využít. První z funkcí (ruleta) generuje kladná celá čísla od 1 do length(cutpoints), každé z nich s pravděpodobností úměrnou hodnotě příslušného prvku vektoru cutpoints. Druhá funkce je pro generování nového bodu heuristikou es-pop. Kromě toho je ještě přidán výsek textu funkce, ve kterém je zapsán výběr heuristik podle četnosti jejich dosavadní úspěšnosti. 58

function [res, 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(1)=cutpoints(1); for i=2:h cp(i)=cp(i-1)+cutpoints(i); end cp=cp/ss; res=1+fix(sum(cp
59

... while (fmax-fmin > my_eps)&(func_evals < d*max_evals) %main loop [hh p_min]=roulete(ni); if p_min<delta cni=cni+ni-n0; ni=zeros(1,h)+n0; nrst=nrst+1; end %reset switch hh % random number of heuristic case 1 y=es_pop(P,0.5); case 2 y=refl_rw(P,2); case 3 y=refl_rw(P,4); case 4 y=refl_fw(P,1); case 5 y=refl_rwd(P,[2 0.5]); case 6 y=refl_rwd(P,[4 0.5]); case 7 y=decrs_rnd(P,[0.5 0.5]); case 8 y=decrs_radp(P,[0.5 0.5 fmin fmax]); end y=zrcad(y,a,b); % perturbation fy=feval(fn_name,y); func_evals=func_evals+1; if fy < fmax % trial point y is good P(indmax,:)= [y fy]; success=success+1; ni(hh)=ni(hh)+1; % zmena prsti qi [fmax, indmax]=max(P(:,d+1)); [fmin, indmin]=min(P(:,d+1)); end end % main loop - end ...

60

Shrnutí

• Soutěž heuristik v algoritmu CRS, pravidla pro posuzování dosavadní úspěšnosti • Soutěž heuristik v algoritmu evolučního prohledávání • Podmínky konvergence algoritmu Kontrolní otázky

1. Co rozumíme pojmem „heuristikaÿ v algoritmech se soutěží heuristik? 2. V čem se liší pravděpodobnosti výběru heuristiky podle vztahu (12.1) a (12.3)? 3. Jakou kombinaci heuristik navrhnete, aby algoritmus splňoval podmínku konvergence? Korespondenční úloha Naprogramujte v Matlabu algoritmus CRS se soutěží nejméně tří vámi zvolených heuristik a porovnejte je alespoň na dvou testovacích funkcích (viz kap. 13) s algoritmy CRS, ve kterých je každá z heuristik užita sama. Pošlete stručnou zprávu o výsledcích testů a m soubory obsahující ověřovaný algoritmus se soutěžícími heuristikami.

61

13

Experimentální ověřování a porovnávání algoritmů

Průvodce studiem V této kapitole je uvedeno několik funkcí, které se užívají v testování stochastických algoritmů pro globální optimalizaci. Z nich si vyberte testovací funkce pro vaše experimenty. Druhá část kapitoly se zabývá testováním stochastických algoritmů. K orientaci v této kapitole budete potřebovat asi dvě až tři hodiny.

13.1

Testovací funkce

K testování stochastických algoritmů pro globální optimalizaci se užívá celá řada funkcí, u kterých je známé správné řešení, tj. globální minimum (většinou lze najít analyticky). Zde uvádíme několik funkcí, které můžeme užít při experimentální testování vámi implementovaných algoritmů. Další testovací funkce naleznete např. v [1, 2, 27, 37] nebo na Internetu. První De Jongova funkce je jednoduchý rotační paraboloid f (x) =

d X

x2i ,

i=1

při experimentech prohledávaný prostor pro d = 3 je obvykle vymezen podmínkou xi ∈ [−5.12, 5.12]. Je to jednomodální funkce, globální minimum je x∗ = (0, 0, 0), f (x∗ ) = 0. Je považována za velmi snadnou úlohu globální optimalizace, pokud algoritmus selhává nebo pomalu konverguje v této úloze, tak nezaslouží další pozornosti. Graf této funkce pro d = 2 je na následujícím obrázku.

62

Druhá De Jongova funkce (Rosenbrockovo sedlo, známá i pod názvem banánové údolí) je definována takto d−1 X   f (x) = 100(x2i − xi+1 )2 + (1 − xi )2 , i=1

prohledávaný prostor je obvykle vymezen podmínkou xi ∈ [−2.048, 2.048], globální minimum je x∗ = (1, 1, . . . 1), f (x∗ ) = 0, obvykle v testech užívané dimenze jsou d = 2, d = 10 a d = 30. Ač je to jednomodální funkce, nalezení minima iteračními algoritmy není jednoduchá úloha, neboť minimum leží v zahnutém údolí s velmi malým spádem a některé algoritmy končí prohledávání předčasně. Graf této funkce pro d = 2 je na následujícím obrázku.

Ackleyho funkce je 

v  u d u1 X x2  − exp f (x) = −20 exp −0.02 t d i=1 i

! d 1X cos 2πxi + 20 + exp(1), d i=1

xi ∈ [−30, 30], x∗ = (0, 0, . . . 0), f (x∗ ) = 0, v testech užívané dimenze d = 2, d = 10 a d = 30. Je to multimodální funkce s mnoha lokálními minimy, středně obtížná. Graf této funkce pro d = 2, xi ∈ [−3, 3], je na následujícím obrázku.

63

Griewangkova funkce,   d d X Y xi xi 2 f (x) = − cos √ + 1, 4000 i=1 i i=1 xi ∈ [−400, 400], x∗ = (0, 0, . . . 0), f (x∗ ) = 0, v testech užívané dimenze d = 2 a d = 10. Nalezení minima této multimodální funkce je považováno za obtížnou úlohu globální optimalizace. Graf této funkce pro d = 2, xi ∈ [0, 15], je na následujícím obrázku.

Rastriginova funkce f (x) = 10d +

d X

[xi 2 − 10 cos(2πxi )]

i=1

xi ∈ [−5.12, 5.12], x∗ = (0, 0, . . . 0), f (x∗ ) = 0, v testech nejčastěji užívaná dimenze je d = 5. Nalezení minima této multimodální funkce je považováno za obtížnou úlohu globální optimalizace. Graf této funkce pro d = 2 je na následujícím obrázku. 64

Schwefelova funkce f (x) =

d X

p xi sin( | xi |),

i=1

xi ∈ [−500, 500], x∗ = (s, s, . . . s), s = 420.97, f (x∗ ) = −418.9829d, v testech užívaná dimenze d = 10. Nalezení minima této multimodální funkce je považováno za středně obtížnou úlohu globální optimalizace. Graf této funkce pro d = 2 je na následujícím obrázku.

Implementace těchto testovacích funkcí je v Matlabu velmi jednoduchá, jak vidíte v nasledujících zdrojových textech algoritmů. 65

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=griewangk(x) % Griewangk’s function, <-600,600>, glob. min f(0)=0 d=length(x); a=sum(x.*x)/4000; j=1:d; b=prod(cos(x./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)); function y=schwefel(x) % Schwefel’s function, <-500,500>, % glob. min f(x_star)=-418.9829*d, x_star=[s,s,...,s], s=420.97 d=length(x); sz=size(x); if sz(1)==1 x=x’; end y=-sum(x.*sin(sqrt(abs(x))));

66

13.2

Experimentální ověřování algoritmů

Numerické experimenty jsou nutné pro ověřování a vzájemné porovnávání stochastických algoritmů. Při porovnávání dvou či více algoritmů je samozřejmě nutné, aby testy byly provedeny za stejných podmínek, tj. zejména při stejně vymezeném prohledávaného prostoru D a při stejné podmínce pro ukončení prohledávání. Základními veličinami pro porovnání algoritmů jsou časová náročnost prohledávání a spolehlivost nalezení globálního minima. Časová náročnost algoritmu se ohodnocuje počtem vyhodnocení účelové funkce v průběhu hledání, takže výsledky jsou srovnatelné bez ohledu na rychlost počítačů užitých k testování. Počet vyhodnocení účelové funkce potřebný k dosažení podmínky ukončení označíme nfe. U testovacích funkcí, kdy řešení problému je známé, se někdy také sleduje počet vyhodnocení účelové funkce nfe1 potřebný k tomu, aby nejlepší bod populace měl hodnotu nižší než je zadaná hodnota fnear , o které se ví, že její dosažení znamená, že algoritmus nejlepším nalezeným bodem se dostatečně přiblížil globálnímu minimu. Někdy se tato hodnota nazývá také VTR, což je zkratka anglického termínu value to reach. Podmínka ukončení bývá formulována jako fmax − fmin < ε. Je však praktické formulovat podmínku ukončení tak, abychom předešli abnormálním ukončení programu např v situaci, kdy prohledávání dojde do nekonečného cyklu. Proto je vhodné přidat další vstupní parametr algoritmu maxevals, což je maximální dovolený počet vyhodnocení funkce na jednu dimenzi prohledáváného prostoru a podmínku ukončení pak formulovat jako fmax − fmin ≤ ε ∨ nfe ≥ maxevals ∗ d, to znamená, že hledání pokračuje tak dlouho, dokud funkční hodnoty v populaci se liší více než požadujeme nebo dokud není dosažen nejvyšší dovolený počet vyhodnocení účelové funkce: fmax − fmin > ε ∧ nfe < maxevals ∗ d. Místo fmax lze v podmínce ukončení užít i jiný kvantil funkční hodnoty v populaci, např. medián. Tím učiníme podmínku ukončení snadněji dosažitelnou, stačí jen, aby např. polovina populace měla funkční hodnoty dostatečně blízké.

67

Z uvedených úvah je zřejmé, že pak musíme rozlišovat následující čtyři typy ukončení prohledávání: • Typ 1 – korektní ukončení, tj. fmax − fmin < ε a fmin ≤ fnear , algoritmus tedy splnil podmínku ukončení před dosažením maximálního dovoleného počtu iterací a současně se dostatečně přiblížil globálnímu minimu. • Typ 2 – pomalá konvergence, fmin ≤ fnear , ale prohledávání bylo ukončeno dosažením maximálního dovoleného počtu iterací. • Typ 3 – předčasné konvergence (premature convergence), fmax −fmin < ε, ale fmin > fnear , tzn. nebylo nalezeno globální minimum, algoritmus buď skončil prohledávání v lokálním minimu nebo se body populace přiblížily k sobě natolik, že jejich funkční hodnoty jsou velmi blízké. • Typ 4 – úplné selhání, prohledávání bylo ukončeno dosažením maximálního dovoleného počtu iterací, aniž byl nalezen bod blízký globálnímu minimu. Po prohledávacích algoritmech globální optimalizace chceme, aby uměly nacházet globální minimum x∗ co nejrychleji, co nejspolehlivěji a aby prohledávání uměly ukončit ve vhodném stadiu prohledávání. Úspěšné prohledávání je tedy jen takové, jehož ukončení je typu 1. Spolehlivost (reliability, R) nalezení globálního minima můžeme charakterizovat jako relativní četnost ukončení typu 1 n1 R= , n kde n1 je počet ukončení typu 1 v n nezávislých opakování. Veličina n1 je náhodná veličina n1 ∼ Bi(n, p) a R je nestranným odhadem parametru p, tj. pravděpodobnosti, že algoritmus v daném problému najde globální minimum (přesněji, že konec prohledávání je typu 1). Rozptyl var(R) = p(1−p) . Pro větší hodnoty n můžeme rozdělení R považovat za přibližně n normální, R ∼ N p, p(1−p) . Pak 100 (1−α)-procentní oboustranný interval n spolehlivosti pro p je [R − 4, R + 4], kde

r α R(1 − R) 4 = tn−1 1 − 2 n V následující tabulce jsou uvedeny hodnoty 4 pro α = 0.05 a pro různá n a R. Pokud byste vyjadřovali reliabilitu R v procentech, pak je potřeba hodnoty 4 v tabulce vynásobit stovkou. Z uvedené tabulky je zřejmé, že pro porovnávání spolehlivosti algoritmů potřebujeme velký počet opakování n. Např. pro n = 100 je možmo za významný považovat až rozdíl v reliabilitě větší než 10 až 20 procent (závisí na hodnotě R). 

68

Hodnoty 4 pro α = 0.05 R n

0.50

0.60

0.70

10 0.358 0.350

0.80

0.90

0.95

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ýsledky testování by tedy vždy měly obsahovat specifikaci testovaných algoritmů včetně hodnot jejich vstupních parametrů, specifikaci testovacích funkcí, počet opakování experimentů, a tabulku s průměrnými hodnotami charakterizující časovou náročnost, vhodné charakteristiky variability a spolehlivost R, která může být uvedena v procentech. Pro názorné porovnání časové náročnosti jsou vhodné krabicové grafy, ve kterých snadno vizuálně porovnáme jak charakteristiky polohy, tak variability. Pokud rozdíly v charakteristickách algoritmů nejsou zjevné z popisné statistiky, použijte vhodné statistické testy. Pro porovnání časové náročnosti to může být dvouvýběrový t-test nebo analýza rozptylu nebom jejich neparametrické analogie. Pokud chcete porovnávat spolehlivosti dvou algoritmů, které se v hodnotách R liší jen málo, takže statistická významnost rozdílu není zřejmá z hodnot pro intervaly spolehlivosti v tabulce, můžete pro čtyřpolní tabulku absolutních četností užít Fisherův exaktní test. Pro některé algoritmy je zajímavé sledovat i jiné veličiny, které charakterizují průběr prohledávacího procesu, např. u algoritmů se soutěžícími heuristikami četnost užití jednotlivých heuristik ap. Volba dalších sledovaných veličin je závislá na testovaném algoritmu. Zajímavou popisnou informací o průběhu vyhledávání je také graf závislosti funkčních hodnot v populaci na počtu vyhodnocení účelové funkce. Průběh této závislosti nám umožňuje posoudit rychlost konvergence algoritmu v různých fázích prohledávání. V takových grafech je ovšem nutné mít charakterisky z více opakování a pokud možno, nejen průměrné hodnoty, ale i minima a maxima.

69

Shrnutí

• Testovací funkce, numerické testování algoritmů • Různé typy ukončení prohledávání • Časová náročnost měřená početem vyhodnocení účelové funkce, spolehlivost nalezení globálního minima • Vyhodnocení experimentu Kontrolní otázky

1. Proč je nutné posuzovat algorimy ne jen z jednoho běhu, ale z více opakování? 2. Jak popsat specifikaci testovacího problému? 3. Jak specifikovat testovaný algoritmus?

70

14

Literatura

[1] Ali M. M., Törn A. (2004). Population set based global optimization algorithms: Some modifications and numerical studies. Computers and Operations Research 31, 1703 – 1725. [2] Bäck, T. (1996). Evolutionary Algorithms in Theory and Practice. Oxford University Press, New York. [3] Bäck, T., Schwefel, H.P. (1993). An Overview of Evolutionary Algorithms for Parameter Optimization, Evolutionary Computation 1, 1–23. [4] Bäck, T., Hoffmeister F., Schwefel, H.-P. (1993). Applications of Evolutionary Algorithms, TR SYS-2/92, University of Dortmund [5] Bäck, T., Fogel, D.B., Michalewicz, Z. (eds). (2000). Evolutionary Computation 2 - Advanced Algorithms and Operators, Institute of Physics Publishing, Bristol and Philadelphia. [6] Dorigo, M., Di Caro, G., Gambardella, L. M. (1999). Ant Algorithms for Discrete Optimization. Artificial Life, 5(2),137–172. [7] Glover F. (1989). Tabu Search - Part I. ORSA Journal of Operational Research 1, 190–206. [8] Glover F. (1990). Tabu Search - Part II. ORSA Journal of Operational Research 2, 4–32 [9] Goldberg, D. E. (1989). Genetic Algorithms in Search, Optimization, and Machine Learning. Reading, Addison Wesley. [10] Holland J. H. (1975). Adaptation in Natural and Artificial Systems, The University of Michigan Press, Ann Arbor, MI. [11] Inberg, L. (1992). Genetic Algorithms and Very Simulated Reannealing: A Comparison, J. Math.Comput. Modeling 16, 87-100. [12] Inberg, L. (1993). Simulated Annealing: Practice versus Theory, J.Math.Comput. Modeling 18, 29-57. [13] Křivý, I., Tvrdík, J. (1995). The Controlled Random Search Algorithm in Optimizing Regression Models. Comput. Statist. and Data Anal., 20(2), 229–234. [14] Křivý, I., Tvrdík, J., Krpec, R. (2000). Stochastic algorithms in nonlinear regression, Comput. Statist. Data Anal. 33, 278–290. [15] Kvasnička, V., Pospíchal, J., Tiňo, P. (2000). Evolučné algoritmy. STU, Bratislava. [16] Lampinen, J. (2001). A Bibliography of Differential Evolution Algorithm. Technical Report. Lappeenranta University of Technology, Department of Information Technology, Laboratory of Information Processing, available in up-to date version via Internet http://www.lut.fi/∼jlampine/debiblio.htm [17] MATLAB, version 6, The MathWorks, Inc., 2000. [18] Michalewicz, Z. (1992). Genetic Algorithms + Data Structures = Evolution Programs. Berlin, Springer Verlag. 71

[19] Michalewicz, Z., Fogel, D.B. (2000). How to Solve It: Modern Heuristics. Springer-Verlag, Berlin Heidelberg New York. [20] Míka, S. (1997). Matematická optimalizace. vydavatelství ZČU, Plzeň. [21] Mišík, L., Tvrdík J., Křivý, I. (2001). On Convergence of a Class of Stochastic Algorithms. In: Antoch, J. and Dohnal, G. (Eds.), Proceedings of ROBUST 2000, JČMF Press, Prague, 198–209. [22] Mühlenbein, H., Schlierkamp-Voosen, D. (1993). The science of breeding and its application to breeder genetic algorithm, Evolutionary Computing 1, 335–360, [23] Nelder, J.A., Mead, R. (1964). A Simplex Method for Function Minimization, Computer J., 7(1) 308–313. [24] Ošmera, P., Kvasnička, V. and Pospíchal, J. (1997). Genetic Algorithms with Diploid Chromosomes, MENDEL’97, 3rd International Conference on Genetic Algorithms, In: MENDEL’97, Technical University, Brno, 111-116. [25] Price, W.L. (1977). A Controlled Random Search Procedure for Global Optimization. Computer J., 20, 367–370. [26] Spall J. C. (2003). Introduction to Stochastic Search and Optimization, Wiley-Intersience. [27] Storn, R., Price, K. (1997). Differential Evolution – a Simple and Efficient Heuristic for Global Optimization. J. Global Optimization, 11, 341–359. [28] Törn, A., Žilinskas A. (1989). Global Optimization, Lecture Notes in Computer Science, No. 350, Springer. [29] Tvrdík J., Křivý I. (1996). A Genetic Type Algorithm for Optimization, MENDEL’96, 2nd International Conference on Genetic Algorithms, In: MENDEL’96, Technical University, Brno, 176-180. [30] Tvrdík, J., Křivý, I., Mišík, L. (2001). Evolutionary Algorithm with Competing Heuristics. In: Proceedings of MENDEL 2001, 7-th Int. Conference on Soft Computing (Matoušek R. and Ošmera P. eds), Technical University Press, Brno, 58–64. [31] Tvrdík, J., Mišík, L., Křivý, I. (2002). Competing Heuristics in Evolutionary Algorithms. In: Sinčák P. et al. (Eds.), 2nd Euro-ISCI Intelligent Technologies - Theory and Applications, IOS Press, Amsterdam, 159– 165. [32] Tvrdík, J., Křivý, I., Mišík, L. (2002). Evolutionary Algorithm with Competing Heuristics in Computational Statistics. In: Härdle, W., Röenz, B. (Eds.), COMPSTAT 2002. Proceedings in Computational Statistics, Physica-Verlag, Heidelberg, 349–354. [33] Wolpert, D. H., Macready, W.G. (1995). No Free Lunch Theorems for Search, SFI-TR-95-02-010, Santa Fe, NM. [34] Wolpert, D. H., Macready, W.G. (1997). No Free Lunch Theorems for Optimization, IEEE Transactions on Evolutionary Computation 1 (1), 67–82. 72

[35] Yang, J.-M., Kao C.-Y. (2000). Integrating adaptive mutations and family competition into genetic algorithms as function optimizer, Soft Computing 4, 89–102. [36] Zelinka, I., Lampinen J. (2000). SOMA – Self-Organizing Migrating Algorithm, In:Proceedings of MENDEL 2000, 6-th Int. Conference on Soft Computing, 177–187. Brno: Technical University. [37] Zelinka, I. (2002). Umělá inteligence v problémech globální optimalizace. BEN, Praha, 2002 WWW stránky: http://math.chtf.stuba.sk/evol/prednaska.htm Stránka k přednáškám prof. Kvasničky a doc. Pospíchala, obsah jejich přednášek je vsouboru sylaby.rtf nebo sylaby.pdf http://alife.santafe.edu:80/ joke/encore/ Nejnovější a dost úplný obsáhlý přehled výsledků, propojení na mnoho dalších stránek o EA a GO, přístup ke zdrojovým textům článků a zpráv http://solon.cma.univie.ac.at/ neum/glopt.html Obsáhlé stránky o GO a EA, mnoho dalších užitečných odkazů na související témata (včetně matematiky a statistiky, adresy stránek mnoha autorů článků a knih věnovaných GO, EA atd.) http://ls11-www.informatik.uni-dortmund.de/people/baeck/ea mj. i postskriptové soubory článků a zpráv, (Baek, Schwefel a spol.) http://www.cs.sandia.gov/opt/survey/main.html přehled, odkazy na základní články http://www.maths.adelaide.edu.au/Applied/llazausk/alife/realfopt.htm několik populárních funkcí pro testování GO, s výsledky (počet vyhodnocení objective function) http://solon.mat.univie.ac.at/ vpk/math/funcs.html 50 testovacích funkcí pro algoritmy GO http://www.lut.fi/ jlampine/debiblio.htm přehled vývoje a aplikací diferenciální evoluce

73