Rozdział 3
Implementacja komputerowa algorytmu genetycznego

Zetknąwszy się po raz pierwszy z algorytmami genetycznymi, wielu użytkowników waha się, nie wiedząc jak i od czego zacząć. Z jednej strony taka niechętna reakcja może dziwić. Przecież, jak widzieliśmy w dwóch pierwszych rozdziałach, algorytmy genetyczne są technicznie całkiem proste, nie wymagają niczego prócz generowania liczb losowych oraz kopiowania ciągów kodowych i wymiany ich fragmentów. Z drugiej strony ta absolutna prostota sama stanowi problem dla wielu użytkowników i programistów działających w sferze biznesu, nauki lub techniki, są to bowiem ludzie nawykli do używania i pisania wysoce zaawansowanych programów, opartych na złożonych podstawach matematycznych, korzystających z powiązanych wzajemnie baz danych i wykonujących zawiłe obliczenia. W dodatku ludzie ci czują się najpewniej, mając do czynienia z budzącą zaufanie powtarzalnością deterministycznych programów komputerowych. Bezpośrednie działania na ciągach bitów, projektowanie specjalnych programów, a nawet losowy charakter operacji wykonywanych przez algorytm genetyczny - wszystko to stanowi serię przeszkód trudnych do przebycia.
Aby pokonać powyższe trudności, skonstruujemy w tym rozdziale struktury danych i algorytmy niezbędne do implementacji opisanego wcześniej elementarnego algorytmu genetycznego. W szczególności, napiszemy w Pascalu program o nazwie SGA (simple genetic algorithm}, realizujący reprodukcję, krzyżowanie i mutację w ramach rozłącznych pokoleń ciągów kodowych - w zastosowaniu do optymalizacji prostej funkcji jednej zmiennej, zakodowanej w układzie dwójkowym. Przedyskutujemy także kilka zagadnień dotyczących implementacji, takich jak dyskretyzacja parametrów, rodzaje kodów, wprowadzanie więzów i przekształcanie funkcji przystosowania, z którymi można się spotkać w rozmaitych zastosowaniach algorytmów genetycznych.
I
76
3. Implementacja komputerowa algorytmu genetycznego
3.1. Strukturydanych
Algorytmy genetyczne działają na populacjach złożonych z ciągów kodowych. Nic wiec dziwnego, że podstawową strukturę danych dla algorytmu genetycznego stanowi populacja ciągów. Istnieją rozmaite sposoby implementacji takich populacji. Dla potrzeb SGA wybierzemy najprostszą metodę: zrealizujemy populację w postaci tablicy osobników [individual}, przy czym każdy osobnik składa się z fenotypu fy>henotype], tj. odkodowa-nego parametru (parametrów), genotypu [genotype], tj. sztucznego chromosomu (ciągu kodowego), wskaźnika przystosowania (wartości funkcji celu) ffiness] oraz pewnych pomocniczych informacji. Na rys. 3.1 przedstawiono schematyczną postać populacji. Fragment programu w Pascalu z wyd. 3.1 zawiera definicję typu reprezentującego populację {populatiori\ w omawianym modelu. Nawet bez formalnej znajomości Pascala, Czytelnik nie powinien mieć większych trudności z uchwyceniem istoty tej specyfikacji.
Numer osobnika Osobniki
Ciąg kodowy X f(X> Inne
1 01111 15 225
2 01001 9 81
3

n 00111 7 49
Rys. 3.1. Struktura populacji ciągów kodowych w algorytmie genetycznym
Na wydruku 3.1 widzimy definicje dwóch stałych: maksymalnego rozmiaru populacji - maxpop i maksymalnej długości ciągu kodowego - maxstring. Narzucają one górne ograniczenia na wielkość populacji i długość ciągu. Po definicjach stałych następują definicje samej populacji i jej składowych (po słowie kluczowym type). Jak widać, typ population stanowi tablicę elementów typu individual (o zakresie od 1 do maxpop). Typ individual stanowi rekord złożony z pól chrom (typu chromosome},fitness (typu real) oraz x (typu real). Reprezentują one odpowiednio sztuczny chromosom, wskaźnik przystosowania i odkodowany parametr". Z kolei typ chromosome sam jest tablicą (o zakresie od 1 do maxstring) elementów typu allele\ ten ostatni jest w tym przypadku synonimem typu boolean (obejmującegojednobitowe wartości true ifalse).
'' Tj. wartość liczbową odpowiadającą chromosomowi (przyp. tlum.).
3.1. Struktury danych
77
const maxpop
maxstring
type
allele
chromosome
individual
- lQO;
- 30; ,. , ....-".'
- boolean; ( Allele - bit position )
- array[l..maxstring] of allele; ( String of bits )
- record
chrora:chroraosome; ( Genotype - bit string ) �
x:real; { Phenotype - unsigned integer )
fitness:real; ( Objective function value }
parentl, parent2, xsite:integer; ( parents & cross pt end; population -array[l..maxpop] ofindividual;
Wyd.3.1. Elementarnyalgorytmgenetyczny,SGA,definicjestrukturdanychwPascalu ,.. r.
W programie SGA operacje genetyczne stosujemy w każdym pokoleniu do całej populacji, jak widać na rys. 3.2. Aby uprościć tworzenie potomków i zastępowanie rodziców, program korzysta z dwóch tablic (oldpop i newpop} odpowiadających rozłącznym populacjom. Ich deklaracje, wraz z deklaracjami innych zmiennych o zasięgu globalnym, podane są w wyd. 3.2. Mając do dyspozycji dwie tablice, możemy bez trudu - używając operacji genetycznych - utworzyć potomstwo elementów oldpop, umieścić je w newpop, a następnie, gdyjużjesteśmy gotowi - skopiować newpop do oldpop. Istnieją inne, bardziej oszczędne pamięciowo metody operowania na populacjach. Moglibyśmy używać jednej tablicy z "mieszaną" populacją, zwracając pilną uwagę na to, kto zastępuje kogo w kolejnych pokoleniach. Nie ma też szczególnego powodu, aby utrzymywać stałą wielkość populacji. Populacje występujące w przyrodzie z całą pewnością zmieniają swą wielkość. Także w przypadku zastosowania algorytmów genetycznych zmiany rozmiaru populacji z pokolenia na pokolenie mogłyby okazać się korzystne. Jednak tutaj bardziej zależy nam na zachowaniu maksymalnej prostoty i dlatego wybraliśmy rozłączne populacje o stałej wielkości. W późniejszych rozdziałach, poświęconych tematyce maszyn uczących się, będziemy musieli jeszcze raz zmierzyć się z problemami dotyczącymi populacji.
Mając zaprojektowane i zaprogramowane struktury danych, powinniśmy teraz zająć się operacjami genetycznymi decydującymi o działaniu SGA - reprodukcją, krzyżowaniem
Pokolenie
T + 1
Reprodukcja
Mutacja
Krzyżowanie
Rys. 3.2. Przejście od starego do nowego pokolenia w SGA
78
3. Implementacja komputerowa algorytmu genetycznego
var oldpop, newpop:population;
( Two non-overlnpping populations )
popsize, lchroro, gen, maxgen:integer; ( Integer global variables )
pcross, pmutatton, sumfitness:real; ( Real global variables )
nmutation, ncross:integer; ( Integer statistics )
avg, max, min:real; ( Real statistics )
Wyd. 3.2. SGA, deklaracje zmiennych globalnych w Pascalu
i mutacją. Przedtem jednak musimy omówić kilka ważniejszych zmiennych o charakterze globalnym wpływających na przebieg wykonania całego programu. Patrząc raz jeszcze na wyd. 3.2, możemy zauważyć pewną liczbę zmiennych typu integer. Są wśród nich zmiennepopsize, lchrom \ gen. Reprezentują one odpowiednio wielkość populacji (którą oznaczaliśmy n), długość ciągu kodowego (/) oraz numer pokolenia (t). Prócz tego, zmienna maxgen określa górną granicę liczby pokoleń. Na wydruku 3.2 widać także kilka zmiennych typu real: pcross, pmutation, sumfitness, avg, max i min. Zmienne pcross i pmutation określają prawdopodobieństwa krzyżowania i mutacji ^)c i pm). Zmienna sumfitness reprezentuje sumę wskaźników przystosowania w populacji CLf)- Odgrywa ona ważną rolę w selekcji proporcjonalnej (wg reguły ruletki). Paru zmiennych globalnych jeszcze nie omówiliśmy; pełny wydruk programu został zamieszczony w dodatku B.
3.2. Reprodukcja, krzyżowanie, mutacja
Każda z trzech operacji prostego trzyczęściowego algorytmu może być zrealizowana w postaci nieskomplikowanej procedury. Nie ma w tym nic dziwnego, wszak głośno zachwalaliśmy prostotę mechanizmu działania tych operacji. Nim zajrzymy w głąb każdej procedury, przypomnijmy sobie wspólną cechę wszystkich trzech operacji: zależność od losowego wyboru. W podanych niżej podprogramach zakłada się istnienie trzech następujących funkcji:
,. .. '.*�: ..��. � . -'j
random - przyjmujejako wartość liczbę pseudolosową z przedziału [0,1] (zmienna losowa o rozkładzie jednostajnym na [0,1]); flip - przyjmuje z zadanym prawdopodobieństwem wartość boolowską true (model
próby Bernoulliego);
rnd - przyjmuje z jednakowym prawdopodobieństwem wartości całkowite między zadanymi liczbami.
Operacja reprodukcji została zaimplementowana w SGA w postaci funkcji select, realizującej wyszukiwanie liniowe na tarczy ruletki wykalibrowanej proporcjonalnie do wskaźników przystosowania ciągów kodowych. W wyniku wywołania funkcja select (wyd. 3.3) przyjmuje wartość indeksu odpowiadającego wybranemu osobnikowi. Zmienna rzeczyv/\slapartsum służy do akumulacji sum częściowych wskaźników przystosowania, wykorzystywanych w tym procesie. Zmienna rzeczywista rand wskazuje położenie, w którym tarcza zatrzymała się po wykonaniu losowego obrotu wokół osi, określonego przez instrukcję
3.2. Reprodukcja, krzyżowanie, mutacja
79
function select(popsize:integer; sumfitness:real;
var pop:population):integer;
( Select a single individual via roulette wheel selection ) var rand, partsum:real; ( Random point on wheel, partial sum }
j:integer; ( population index ) begin
partsum :- 0.0; j :- 0; ( Zero out counter and accumulator ) rand : random * sumfitness; ( Wheel point calc. uses random number [0,1] ) repeat { Find wheel slot ) j :- j + 1;
partsum : partsum + pop[j].fitness; i
until (partsum > rand) or (j popsize); { Return individual number ) select :- j; end;
Wyd. 3.3. Funkcja paskalowa se/ecfimplementuje selekcję wg reguły ruletki
rand := random * sumfitness ,
Suma wszystkich wskaźników przystosowania (obliczona w procedurze statistics) zostaje tu pomnożona przez znormalizowaną liczbę pseudolosową wygenerowaną przez random. Następnie w pętli repeat-until poszukuje się pierwszej sumy częściowej nie mniejszej od krytycznej wartości rand. Wywołanie kończy się przypisaniem funkcji select wartości bieżącego indeksu osobnika/
Jest to zapewne najprostsza metoda implementacji wyboru z populacji. Istnieją bardziej efektywne sposoby urzeczywistnienia tej operacji (wyszukiwanie binarne z pewnością przyspieszy bieg spraw), a także wiele innych metod tendencyjnego wyboru zakładających uprzywilejowanie najlepszych osobników. Z niektórymi z nich zapoznamy się w następnych rozdziałach, na razie jednak poprzestaniemy na tym podstawowym mechanizmie.
Podprogram select dostarcza nam prostego narzędzia wyboru osobników do następnego pokolenia'\ Jak pamiętamy z wcześniejszego opisu, następnym krokiem jest krzyżowanie. W SGA operacja krzyżowania została zrealizowana w postaci procedury o nazwie crossover (wyd. 3.4). Podprogram ten pobiera dwa rodzicielskie ciągi kodowe parentl oraz parent2 i wytwarza dwa ciągi potomne childl oraz child2. Do procedury są przekazywane prawdopodobieństwa krzyżowania pcross i mutacji pmutation, a także długość ciągu kodowego lchrom, licznik krzyżowań ncross i licznik mutacji nmutation.
Instrukcje wewnątrz procedury crossover odzwierciedlają opis operacji krzyżowania podany w rozdziale 1. Na wstępie podejmujemy decyzję, czy będziemy dokonywać skrzyżowania danej pary chromosomów rodzicielskich. Konkretniej, wykonujemy rzut tendencyjną monetą, dla której orzeł (true) wypada z prawdopodobieństwem pcross. Doświadczenie tojest symulowane za pomocą wywołania funkcji boolowskiej/%, która z kolei wywołuje generator liczb pseudolosowych random. Jeżeli decyzja jest pozytywna, trzeba wybrać punkt krzyżowania między 1 a ostatnim możliwym miejscem. Punkt krzyżowania zostaje wylosowany za pomocą wywołania funkcji rnd, która generuje
" Właściwie do puli rodzicielskiej Q)rz.yp. ttum.).
80
3. Implementacja komputerowa algorytmu genetycznego
pseudolosową liczbę całkowitą zawartą między dwiema wskazanymi liczbami (tj. 1 i lchrom - 1). Jeżeli decyzja w sprawie krzyżowaniajest negatywna, punkt krzyżowania otrzymuje wartość lchrom (czyli leży za chromosomem), dzięki czemu równoległa operacja mutacji zostanie wykonana mimo braku krzyżowania. W końcu, w dwóch pętlach for-do stanowiących zakończenie procedury, dokonuje się częściowa wymiana genów związana z operacją krzyżowania. W pierwszej pętli część materiału genetycznego zostaje skopiowana zparentJ do childl oraz zparent2 do child2. W drugiej pętli pozostały materiał genetyczny zostaje skopiowany z parentl do child2 oraz z parent2 do childl. W obu przypadkach kopiowane geny mogą ulegać mutacji realizowanej za pomocą wywołania boolowskiej (czy też ,,allelicznej") funkcji mutation.
{ Do crossover with p(cross) )
( Cross beCween 1 and 1-1 )
( lncrement crossover counter ) ( Otherwise set cross site to force mutation
procedurę crossover(var parentl, parent2, childl, child2:chromosome; var lchrom, ncross, nmutation, jcross:integer; var pcross, pmutation:real);
( Cross 2 parent strIngs, place in 2 child strings ) var j:integer; begin
if flip(pcross) then begin jcross :- rnd(l,lchrom-l) ncross ;- ncross + 1; end else '�.��',:� jcross :- lchrom;
( lst exchange, 1 to 1 and 2 to 2 ) for j : 1 to jcross do begin
childl[j] : mutation(parentl[j], pmutation, nmutation); .,,-. child2[J] : mutation(parent2[j], pmutation, nmutation); end;
( 2nd exchange, 1 to 2 and 2 to 1 ]
if jcrossOlchrom then ( Skip if cross site is lchrom--no crossover ) for j :- jcross+1 to lchrom do begin
childl[j] :- mutation(parent2[j], pmutation, nmutation); child2[j] : mutation(parentl[j], pmutation, nmutation); end; end;
Wyd. 3.4. Procedura crossoverimplementuje krzyżowanie proste Qednopunktowe)
Funkcja mutation (wyd. 3.5) decyduje o mutacji poszczególnych genów. Korzysta ona z funkcji flip (rzut tendencyjną monetą) w celu określenia, czy należy zmienić true nafalse (1 na 0) lub odwrotnie - czy też nie należy. Oczywiście^fyj - w rezultacie wywołania generatora liczb pseudolosowych random - da orła (tj. wartość true) przecięt-
functiop mutat:ion(alleleval:allele; pmutation:real;
var nmutation:integer):allele;
( Mutate an allele w/ pmutation, count number of mutations ) var mutate:boolean; ; begin
mutate :- flip(pmutation); ( Flip the biased coin ) if mutate then begin
nrautation :- nmutation + 1;
mutation :- not alleleval; ( Change bit value ) end else
mutation :- alleleval; ( No change ) end;
Wyd. 3.5. Funkcja paskalowa mutation implementuje mutację punktową (jednobitową)
3.3. Czas reprodukcji, czas krzyżowania .
81
nie jedynie w pmutation � 100 procentach przypadków. Podprogram rejestruje także liczbę dokonanych mutacji, zwiększając za każdym razem wartość licznika nmutation o 1. Podobnie jak w przypadku reprodukcji jest możliwe ulepszenie implementacji. Na przykład wielu wywołań generatora liczb losowych można by uniknąć, decydując, kiedy ma nastąpić kolejna mutacja, zamiast wywoływać za każdym razem funkcje///p. I znow,jak poprzednio w tym rozdziale, zrezygnujemy z zawiłych udoskonaleń, zadowalając się podstawowym mechanizmem.
I tak oto trzy główne elementy naszej układanki algorytmicznej okazały się nie takie znów trudne do złożenia. W punkcie tym pokazaliśmy, że każdy z nich można łatwo zaprogramować i zrozumieć. W następnym punkcie będziemy kontynuować składanie jeszcze większego "obrazka" algorytmu genetycznego, koordynując wzajemnie reprodukcję, krzyżowanie i mutację w ramach jednego pokolenia.
3.3. Czas reprodukcji, czas krzyżowania
Z taką Wielką Trójką - zaprojektowaną i zaprogramowaną - utworzenie nowej populacji ze starej nie przedstawia już wielkiej trudności. Odpowiedni ciąg instrukcji jest podany w wyd. 3.6, w procedurze generation. Zaczynając od indeksuy'=l i kontynuując po-
procedure generation;
( Create a new generation through select, crossover, and mutation )
( Note: generation assumes an even-numbered popsize )
var j, matel, mate2, jcross:integer;
begin
j :- 1;
repeat ( select, crossover, and mutation until newpop is filled ) matel :- select(popsize, sumfitness, oldpop); ( pick pair of mates } mate2 :- select(popsize, sumfitness, oldpop);
{ Crossover and mutation - mutation embedded within crossover ) crossover(oldpop[matel].chrom, oldpop[mate2].chrom, newpop[j ].chrom, newpop[j + l].chrom, lchrom, ncross, nmutation, jcross, pcross, pmutation);
( Decode string, evaluate fitness, & record parentage data on both children ) with newpop[j ] do begin x :- decode(chrom, lchrom); fitness :- objfunc(x); parentl :- matel; parent2 :- mate2; xsite :- jcross; end;
with newpop[j+l] do begin x :- decode(chrom, lchrom); fitness : objfunc(x); parentl' : matel; ' '-� parent2 :- raate2; xsite :- jcross; end; { Increment population index ) , , �
J :- J + 2;
until j>popsize V \
end;
Wyd.3.6.Procedurageneraf/ontworzynowapopulacjezpoprzedniej ;
82
3. Implementacja komputerowa algorytmu genetycznego
stepowanie aż do przekroczenia wielkości populacji popsize, wybieramy parę partnerów matel \ mate2 za pomocą dwukrotnego wywołania funkcji select. Następnie przeprowadzamy krzyżowanie chromosomów połączone z mutacją, używając do tego procedury crossover (która z kolei zawiera niezbędne wywołania funkcji mutation}. W ostatnim porywie dzikiej aktywności dekodujemy parę chromosomów, obliczamy wskaźniki przystosowania (wartości funkcji celu) i zwiększamy indeks osobnikówy o 2.
Zapoznawszy się szczegółowo z podprogramami select, crossover i mutation, możemy teraz skupić się na dwóch dalszych, zasygnalizowanych wyżej. Są to podprogramy związane z zastosowaniami. Dla każdego zadania musimy napisać procedurę dekodującą ciągi (chromosomy), służącą do odtwarzania wartości parametru (czy też parametrów) charakteryzującego (charakteryzujących) to zadanie. Musimy mieć także procedurę, która wyznacza na podstawie wartości tego parametru (parametrów) wartość funkcji celu związanej z danym zbiorem parametrów. Te dwie procedury, o nazwach decode i objfunc, to jedyne miejsca, w których koła algorytmu genetycznego dotykają szosy zastosowań. Dla różnych zastosowań będą niekiedy potrzebne różne procedury dekodu-jące (chociaż pod koniec tego punktu omówimy pewne standardowe procedury, które okazały się użyteczne w szerszym zakresie), a w każdym razie - różne procedury ewalu-acji funkcji przystosowania. Powiedziawszy to wszystko, warto by jednak skupić się na konkretnej procedurze dekodującej i konkretnej funkcji przystosowania. Aby zachować ciągłość z wcześniejszym materiałem, pozostaniemy przy użyciu zwykłego kodu dwójkowego oraz prostej funkcji potęgowej jako funkcji przystosowania, zwiększając jednak wartość wykładnika do 10 tf(x)=x]n).
function decode(chrom:chromosome; lbits:integer):real;
( Decode string as unsigned binary integer - true-1, false-0 )
var j:integer; :,
accum, powerof2:real; begin
accum :- 0.0; powerof2 :- 1; for j :- 1 to lbits do begin
if chrom[j] then accum :- accum + powerof2; powerof2 :- powerof2 * 2;
end; - ... -. / * , �, : - f ,. . , ...,,..
decode :- accum; end;
Wyd. 3.7. Funkcja paskalowa decode dekoduje ciąg dwójkowy do postaci liczbowej
W programie SGA używa się podprogramu dekodującego zrealizowanego w postaci funkcji paskalowej decode (wyd. 3.7). Chromosom jest dekodowany bit po bicie począwszy od najmłodszego (pozycja 1), przez sumowanie kolejnych potęg dwójki (zmienna poweroftwo} odpowiadających pozycjom zawierającym wartość true. Wywołanie kończy się przypisaniem funkcji decode końcowej sumy, zakumulowanej w zmiennej accum.
Funkcja celu użyta w SGA to prosta funkcja potęgowa, podobna do tej, z którą mieliśmy do czynienia w rozdziale 1. W programie obliczamy wartość funkcji f(x) = (x/coeff)l(\ Wartość współczynnika coe^została wybrana tak, aby znormalizować parametr x w przypadku ciągu kodowego o długości lchrom = 30. Mamy więc
3.4. Program główny
83
functlon objfunc(x:real):real;
( Fitness functian - f(x) - x**n )
const coef - 1073741823.0; ( Coefficient to normalize domaln }
n 10; ( Power of x )
begin objfunc : power( x/coef, n ) end;
Wyd. 3.8. Funkcja paskalowa objfunc oblicza wartość funkcji przystosowania f(x) = cx10 dla zdekodo-wanego parametru x
coeff=2V}- 1 = 1073741823. Po znormalizowaniu wartości x maksimum funkcji/(*) wypada w punkcie x=2y)\ i wynosi 1 (dla lchrom = 30). Na wydruku 3.8 jest podana przejrzysta implementacja funkcji celu w postaci funkcji paskalowej objfunc.
3.4. Program główny
Zaprojektowaliśmy struktury danych. Zaprogramowaliśmy operacje genetyczne. Umie-rny dekodować ciągi i obliczać wskaźniki przystosowania. Teraz nadszedł czas, aby okleić ten pakunek taśmą, sprawdzić wytrzymałość i przekazać do dalszego użytku. Wydruk 3.9 przedstawia program główny dla SGA. Zaczyna się on dość niewinnie, od ustawienia licznika pokoleń na 0, gen := 0. Stopniowo dodajemy pary, wczytując parametry programu, generując losową populację, obliczając statystyki dla populacji początkowej i drukując specjalny raport początkowy - za pomocą procedury initialize. Nie będziemy się tu zatrzymywać nad szczegółami inicjalizacji. Zainteresowany Czytelnik może odwołać się do dodatku B, gdzie jest zamieszczona pełna treść programu SGA.
begin ( MaIn program ) gen :- 0; l Set things up ) initialize; repeat ( Main iterative loop ) p
gen :- gen + 1;
generation;
statistics(popsize, max, avg, min, sumfitness, newpop);
report(gen);
oldpop :- newpop; ( advance the generation ) until (gen > maxgen) end. ( End main program )
Wyd. 3.9. Program główny w SGA
Po długich przygotowaniach docieramy wreszcie do głównej pętli programu (instrukcja repeat-until). Teraz raz za razem zwiększamy licznik pokoleń, tworzymy nowe pokolenie za pomocą procedury generation, obliczamy aktualne statystyki w procedurze statistics, drukujemy raport o bieżącym pokoleniu za pomocą procedury report i jednym zamachem dokonujemy wymiany pokolenia:
; .��. �. < ,ł , . � *
oldpop :=newpop;
Proces ten jest kontynuowany nieustająco, krok po kroku, do chwili, gdy licznik pokoleń przekroczy maksimum, zmuszając maszynerię do gwałtownego zatrzymania.
84
3. Implementacja komputerowa algorytmu genetycznego
procedurę statistics(popsize:integer;
var raax,avg,rain,sumfitness:real;
var pop:population); ( Calculate population statistics ) var j:integer; begin ( Initialize )
sumfitness rain
- pop[l].fitness;
- pop[l].fitness;
- pop[1J.fitness;
( Loop for raax, min, sumfitness )
for j :- 2 to popsize do with pop[j] do begin
( Accumulate fitness sum ) ( New max ) ( New min )
sumfitness :- sumfitness + fitness; if fitness>max then max :- fitness; if fitness( Calculate average ) : avg : sumfŁtness/popsize; end;
Wyd. 3.10. Procedurasfafef/ćswyznaczaważne statystyki populacyjne
W pogoni za ogólnym obrazem zapomnieliśmy o kilku istotnych szczegółach. Podprogram statystyczny statistics (wyd. 3.10) oblicza średni, maksymalny i minimalny wskaźnik przystosowania, a także wartość zmiennej sumfitness potrzebną w procedurze selekcji. Przedstawiona wersja procedury statistics może być uznana za minimalne akceptowalne rozwiązanie. Moglibyśmy (i zapewne powinniśmy) śledzić wiele innych interesujących statystyk populacyjnych. Często, na przykład, praktykuje się zbieranie statystyk dotyczących zbieżności poszczególnych alleli. Można także rejestrować najlepszy zaobserwowany ciąg kodowy albo k najlepszych ciągów. Do bardziej szczegółowych analiz mogą być przydatne odchylenia standardowe lub nawet histogramy populacji. Dzięki zgrupowaniu wszystkich funkcji statystycznych w obrębie jednej procedury można łatwo uzupełnić ją o dowolne z wymienionych obliczeń.
Procedura report (wyd. 3.11) drukuje pełny raport na temat populacji, włączając w to ciągi kodowe, wskaźniki przystosowania i wartości parametru x. I znowu nie ulega wątpliwości, że w pracy z algorytmem genetycznym przydałby się bogaty zestaw opcji do tworzenia tabel i wykresów. Prosta procedura report jest wygodnym narzędziem, gdyż umożliwia porównanie - wiersz po wierszu - następujących po sobie pokoleń. To z kolei daje możliwość sprawdzenia operacji genetycznych i analizy zdarzeń prowadzących do konstrukcji najlepszych osobników.
3.5. Czy wszystko działa jak należy?
Odbyliśmy mozolny marsz, krok po kroku, centymetr po centymetrze, przez instrukcje programu SGA, zdobywając w ten sposób lepsze pojęcie o tajnikach programowania algorytmu genetycznego. Oczywiście, niejedna droga prowadzi do celu; istnieją publicznie dostępne i działające efektywnie programy optymalizacyjne dla różnych dziedzin problemowych, wyposażone w rozliczne "wodotryski" (Booker i De Jong, 1985; De Jong, 1982; Grefenstette, 1984a, 1984b). My również będziemy stopniowo dodawali różne ważne
3.5. Czy wszystko działa jak należy?
85
{ report.sga: contains writechrom, report )
procedurę writechrom(var out:text; chrom:chromosome; lchrora:integer); ( WriCe a chromosome1 as a string of l's (true's) and O's (false's) ) var j:integer; begin for j :- lchrom downto 1 do
if chrom[j] then write(ouC,'l') else write(out,'0'); end;
procedurę report(gen:integer); ( Write the population report ) const linelength 132; var j:integer; begin
repchar(lst,
repchar(lst,
repchar(lst,
repchar(lst,
writeln(lst)
#
,linelength); writeln(lst); ,50); writeln(lst,'Population Report'); ,23); write(lst,'Generation ',gen-l:2); ,57); writeln(lst,'Generation ',gen:2);
write(lst, # string x fitness');
write(lst, # parents xsite');
writeln(lst, ' string x fitness');
repchar(lst,'-',linelength); writeln(lst); for j : 1 to popsize do begin write(lst,j:2, ') '); ( 01d string ) ' ~
with oldpop[j] do begin
writechrom(lst,chrom,lchrom);
write(lst,' ', x:10, ' ', fitness:6:4, � |'); ,-,\: ^;
end;
( New string ) with newpop[j] do begin
write(lst,' ', j:2, ') (', parentl:2, ',', parent2:2, ')
xsite:2,' '); writechrom(lst,chrora,lchrom); writeln(lst, ' ',x:10,' ', fitness:6:4); end;
end; "' ;�'� i ' �
repchar(lst,'-',linelength); writeln(lst); ( Generation statistics and accumulated values ) writeln(lst,' Note: Generation ', gen:2, ' & Accumulated Statistics: '
,' max', max:6:4,', min-', min:6:4, ', avg*', avg:6:4, ', sum-' ,sumfitness:6:4, ', nmutation-', nrautation, ', ncross- ', ncross); repchar(lst,'-',linelength); writeln(lst); page(lst); end; ' !
Wyd. 3.11. Procedury report i writechrom przygotowują raporty o stanie populacji
uzupełnienia, w tym i w następnych rozdziałach. Odrzućmy więc pokusę i powstrzymajmy się na razie od wprowadzania fantazyjnych ulepszeń. W tym punkcie chcemy się zająć samym szkieletem algorytmu genetycznego i zbadać, czy działa on jak należy.
Określiliśmy już proste zadanie kontrolne. Ciąg kodowy reprezentuje 30-bitową liczbę całkowitą. Funkcja przystosowania/ma postać funkcji potegowej/(*) = (*/c)", gdzie c jest stałą normalizującą x, a jako n wybraliśmy 10. Niektórzy Czytelnicy mogą zaprotestować, dziwiąc się dlaczego zmieniliśmy funkcję przystosowania (w poprzednich rozdziałach rozpatrywaliśmy funkcję/(*)=*2). W rzeczy samej, zmiana ta ma
86
3. Implementacja komputerowa algorytmu genetycznego
utrudnić zadanie algorytmu genetycznego, jak pokazuje rys. 3.3. Przy większym wykładniku średnia wartość funkcji jest mniejsza, a zatem mniejsza cześć dziedziny funkcji odpowiada wartościom większym od pewnej ustalonej wielkości. W konsekwencji losowa populacja początkowa nie będzie zawierać zbyt dobrych punktów startowych; jest to więc mocniejszy test sprawności algorytmu.
f(x)
Rys.3.3.Wykresyfunkcjix2ix10naodcinkujednostkowym
Aby określić precyzyjnie warunki symulacji, wybierzmy eksperymentalne wartości parametrów programu. Studium De Jonga (1975) poświęcone optymalizacji funkcji przy użyciu algorytmów genetycznych sugeruje, na podstawie szeregu badań parametrycznych wykonanych dla zestawu pięciu funkcji różnych typów, że w celu zapewnienia wysokiej sprawności algorytmu genetycznego należy dobierać duże prawdopodobieństwo krzyżowania, małe prawdopodobieństwo mutacji (odwrotnie proporcjonalne do wielkości populacji) oraz umiarkowaną wielkość populacji. Idąc za tymi wskazówkami ustalamy następujące wartości parametrów w pierwszych symulacjach komputerowych:
pmutation = 0,0333 (prawdopodobieństwo mutacji)
pcross = 0,6 (prawdopodobieństwo krzyżowania)
popsize =30(wielkoscpopulacji,n) ..,.,, ...
W rozdziale 1 do odręcznej symulacji algorytmu genetycznego wybraliśmy krótkie ciągi kodowe (krótkie, jak na standardy przyjęte dla algorytmów genetycznych) o długości / = 5. Wynikła stąd śmiesznie mała przestrzeń poszukiwań licząca tylko 25 = 32 punkty, dla której zastosowanie algorytmu genetycznego praktycznie mija się z celem. Jakakolwiek metoda przeglądu czy poszukiwania losowego szybko znalazłaby dobre rozwiązania. Oczywiście chodziło nam wówczas o osiągnięcie pedagogicznej przejrzystości, a rozmiar przestrzeni poszukiwań nie był przedmiotem zainteresowania. Obecnie jednak, kiedy zamierzamy wykonać ostrzejszy test sprawności algorytmu, długość ciągu została zwiększona, podobnie jak wykładnik funkcji przystosowania. Przy długości ciągu
3.5. Czy wszystko działa jak należy?
87
lchrotn = 3Q mamy do czynienia ze znacznie większą przestrzenią poszukiwań i ani błądzenie przypadkowe, ani przeglądanie nie okazałyby sięjuż tak korzystne. Dla lchrom = 3Q mamy 230= 1,07- 10'J punktów. Przy ponad miliardzie punktów w przestrzeni rozwiązań żadna z metod poszukiwania punkt-po-punkcie nie wydaje się mieć szans na szybkie osiągnięcie dobrego wyniku. Ponadto, zwiększenie wykładnika powoduje, że tylko 1,05% punktów osiąga wskaźnik przystosowania powyżej 0,9 (por. rys. 3.3). Wprowadzenie obydwu modyfikacji umożliwi zatem lepsze przetestowanie sprawności algorytmu.
SGA Parameters
Population size (popsize) Chromosome length (lchrom) Maximum # of generation (maxgen) Crossover probability (pcross) Hutation probability (pnutation)
30
30
10
6.0000000000E-01
3.3300000000E-02
Initial Generation Statistics
Initiat population maximum fitness * 2.8241322532E-01
Initial poputation average fitness = 3.4715832788E-02
Initial population minimum fitness � 1.U06151375E-10
Initial population sum of fitness - 1.04U749837EtOO
Wyd. 3.12. Raport początkowy z przykładowego przebiegu programu SGA
Uruchamiamy program i pozwalamy mu działać przez siedem pokoleń. Raport inicjalny dla przebiegu jest podany w wyd. 3.12, a wyd. 3.13 ukazuje pokolenie począt-
1,0
0,5-
o,o
max. avg.
3 4 5
Numer pokolenia
Rys. 3.4. Test kontrolny SGA. Maksymalne (max) i średnie (avg) wskaźniki przystosowania ciągów kodowych w kolejnych pokoleniach
88

# Population Report Generation 0 � ' ' string x fitness # ?'^! ';<;f'; ��.-, Generation 1 parents xsite string x fitness
1) 1 1100001 1001 100000101 1 1 1 1 1 0100 9.4621EOS 0.2824 1) ( 1, 19) 30 111000011001100000101101110100 9.4621E*08 0.2824
2) 110011001011010111000000100001 B.5862E*08 0.1069 2) < 1. 19) 30 110101011110001110010011000101 8 9712E*08 0.1658
3) 010101111001011110000010001101 3.6739E+08 0.0000 3) (23, 19) 19 110101011100000011010110000001 8 9655E*08 0.1647
4) 011001111000011101101111011010 4.3423E*08 0.0001 4) (23, 19) 19 11111 100100000001001001 1000101 1 0591E->09 0.8715
5) 011111111010010101011010110110 5.3539E*08 0.0009 5) (19, 1) 11 111000011001100000110011000101 9 46J1E*08 0.2824
ć) 101101111001000011000101101101 7.6993E*08 0.0359 6) (19, 1) 11 110101011110000010001111110100 8 9707E*08 0.1657
7) 000110101111100001001011111000 1.1312E*08 0.0000 7) (16, 1) 6 111000011001100000101111100100 9 4621E*08 0.2824
8) 010100111010111010001111100010 3.5099E*08 0.0000 8) (16, 1) 6 110011010101101100011001110100 8 6132E*08 0.1103
9) 011011001110001010011110001011 4.5670E+08 0.0002 9) (23, 17) 30 101111001000010011110110010001 7 9071E*08 0.04Ć9
10) 010011010110101000001101011011 3.2470E*08 0.0000 10) (23, 17) 30 110011101001000100011101100011 8 6640E+08 0.1170
11) 010111011010101011001101000010 3.9287E+08 0.0000 11) ( 6, 26) 30 101101110001000011000101101101 7 6783E*08 0.0350
12) 010011010011001010110010001110 3.2379E+08 0.0000 12) ( 6, 26) 30 110010000101010100110110011110 8 4026E08 0.0861
13) 110000100101101110110000100111 8.1520E+08 0.0636 13) ( 2, 19) 10 110101011110100010010000100001 8 9720E+08 0.1659
H) 010110001001111101000100110001 3.7171E+08 0.0000 14) ( 2, 19) 10 110011011011010111000011000101 8 6281E*08 0.1122
15) 010110111110000101101110011010 3.8538E+08 0.0000 15) (17, 17) 26 110011111001000100011100100011 8 7060E*08 0.1228
16) 110011010101100100011001100000 8.6129E->08 0.1103 16) (17, 17) 26 110010111101000100011100100011 8 5487E*08 0.1023
17) 110011111101000100011100100011 8.7165E+08 0.1243 17) (19, 17) 30 1101010111 1000001001001 10001 1 1 8 9707E*OB 0.1657
18) 100000000100010111100101011100 5 ,3802E+08 0.0010 18) (19, 17) 30 110011111101000000011101000011 8 7163E*08 0.1243
19) 1101010111 1000001001001 1000101 8.9707E+08 0.1657 19) (19, 19) 24 110101011110000010010011000101 8 9707E*08 0.1657
20) 010011111000010001000011011101 3.3352E*08 0.0000 20) (19, 19) 24 110101011110000010010011000101 8 9707E+OB 0.1657
21) 001101011100110111010010000011 2.2567E*08 0.0000 21) (26, 17) 30 110010000101010100110111011110 8 4026E+08 0.0861
22) 000110011111000001100100110110 1.0880E+08 0.0000 22) (26, 17) 30 110000111101000100011100100011 8 2132E*08 0.0686
23) 101111001000000011110110010001 7.9064E+08 0.0469 23) (23, 1) 3 111000011001100000001111110001 9 4621E*08 0.2824
24) 011101100001100010100101100111 A.9533E+08 0.0004 24) (23, 1) 3 101111001000000011110100010100 7 9064E*08 0.0469
25) 010110111010001101010001010010 3.8436E*08 0.0000 25) (27 1) 10 110000011001100000101001110011 8.1199E+08 0.0612
26) 110010000101010100110110011110 8.4026E*08 0.0861 26) (27, 1) 10 101001100001010101001101110100 6 9660E+08 0.0132
27) 101001100101110101001001100011 6.9778E+08 0.0134 27) ( 1, 17) 25 110010001001100000101111110100 8 4135E*08 0.0873
28) 100011110111010000100000110010 6.0169Et08 0..0031 28) ( 1, 17) 25 111001111101000100011100100011 9 .7231E08 0.3707
29) 010010100010000100001101100011 3.1092E+08 0.0000 29) ( 1 19) 23 110101010001100000101111110100 8 .9378E*08 0.1597
30) 001011001111001110010101100011 1.8854E+08 0.0000 30) ( 1, 19) 23 111000011110000000010011000101 9 .4739E*08 0.2859
Notę: Generation 1 S Accumulated Statistics: max*0.8715, nin*0. 0132, �vg =0. 1732, sum=5.1967, nmutation=35, ncross= 10
Wyd. 3.13. SGA, raport o pokoleniach 0 i 1
s
kowe (gen = 0) i pierwsze pokolenie potomne zestawione obok siebie. Sredni wskaźnik przystosowania dla pokolenia początkowego wyniósł 0,0347. Nietrudno obliczyć, że taki sam średni wskaźnik dla całego przedziału (przestrzeni rozwiązań) wynosi 0,0909. W pewnym sensie mieliśmy więc pecha (choć nie przesadnie wielkiego) z wyborem populacji początkowej. Biorąc ponadto pod uwagę, że najlepszy wskaźnik przystosowania wyniósł /mui=0,2824, moglibyśmy oczekiwać średnio 30(l-0,2824"'')-3,56, tj. w zaokrągleniu czterech ciągów w 30-elementowej losowej populacji o wyższym wskaźniku przystosowania. Mieliśmy więc nie tylko ogólnego pecha, mieliśmy także pecha z maksimum. Mimo tych niepomyślnych początków, algorytm szybko dochodzi do dobrych wyników, co wyraźnie widać po pierwszej turze reprodukcji, krzyżowania i mutacji. W pierwszym pokoleniu potomnym znajdujemy bardzo dobry ciąg o wskaźniku przystosowania 0,8715. W kolejnych turach następuje dalsza poprawa wyników zarówno pod względem maksymalnego, jak i średniego przystosowania populacji, co widać na rys. 3.4. Pod koniec przebiegu obserwujemy coś w rodzaju zbieżności (wyd. 3.14). Gdybyśmy prze-
89
! . -� ' � i ' ! 1 śledzili zawartości poszczególnych pozycji ciągów kodowych należących do siódmego pokolenia, moglibyśmy zauważyć w większości przypadków znaczną zgodność. Zjawisko to wystąpiło, mimo że nie znaleźliśmy optymalnego rozwiązania; tym niemniej zbliżyliśmy się do niego. W szóstym pokoleniu wystąpił osobnik o wskaźniku przystosowania 0,9807. Jest to wynik bliski optymalnemu, choć nie optymalny (wspomniany punkt zalicza, się do czołowych 0,19% punktów całej przestrzeni). Zbieżność, ale bez gwarancji optymalności niepokoi wiele osób, które podchodzą do algorytmów genetycznych z doświadczeniem wyniesionym ze stosowania innych, tradycyjnych metod optymalizacji.
4 Istnieją sposoby, aby spowolnić tę tak zwaną przedwczesną zbieżność [premature con-
T i vergence}\ niektóre z nich omówimy w tym i następnych rozdziałach. Pozostaje jednak
ł faktem, że algorytmy genetyczne nie dają żadnej gwarancji zbieżności dla zadań dowol-
nego typu. Potrafią one szybko odszukać interesujące obszary przestrzeni rozwiązań, ale
i t należą do metod słabych, pozbawionych gwarancji charakterystycznych dla algorytmów
t i 4 zbieżnych. Nie zmniejsza to jednak ich użyteczności. Przeciwnie, metody gwarantujące
j i Population Report Generation 6 . , Generation 7
ł
ł ' * J string x fitness # parents xsite string x fitness
1> 111100011010010010101111110001 1.0135E+09 0.5615 1) ( 8, 6) 7 111111001011000000101111000100 1 .0599E+09 0.8779
4 2) 111111001000100010011111110011 1.0592E+09 0.8726 2) ( 8, 6) 7 111111001000000010011011101100 1.0591E*09 0.8715
3) 111111001000000010011111110100 .0591E+09 0.8715 3) ( 9,24) 9 111111001000000010011010100111 1.0591E*09 0.8715
ł 4) 111110011110000000101111100100 .0481E+09 0.7849 4) ( 9,24) 9 111111101000000010000111110111 1.0675E+09 0.9430
t 5) 111111001101100001101111010110 .0605E+09 0.8834 5) ( 6,18) 3 111111011000001110101111100100 1.0633E*09 0.9070
6, 111111001011000000101111101100 .0599E+09 0.8779 6) ( 6,18) 3 111110001011000000111111101100 1.0431E*09 0.7484
J 7) 111111001001100000101111101100 .0595Et09 0.8747 7) (10,22) 22 101110011100000010011111110111 7.7910E*08 0.0405
4 8) 111111001000000010011011000100 .0591E+09 0.8715 8) (10,22) 22 111111001111010000101111100100 1.0610E*09 0.8872
9) 111111101000000010000011100111 .0675E+09 0.9430 9) (15,15) 30 111110011010000010011111110001 1.0470E*09 0.7772
] 10) 111111001100000010011111110111 1.0601E+09 0.8801 10) (15,15) 30 111110011010000010011111111001 1.0470E+09 0.7772
1" 111111001000000010011111110111 1.0591E+09 0.8715 11) (12,12) 5 111111111000000000011111010001 1.0716E+09 0.9807
12) 111111101000000010011111110001 1.0675E+09 0.9430 12) (12,12) 5 111111101000000010011111110001 1.0675E+09 0.9430
l > 111011011000000010000011000101 9.9616E+08 0.4724 13) (26,16) 15 111111001010000000101101110111 1.0596E+09 0.8757
14) 111111001001110000101111110000 1.0595E+09 0.8752 14) (26,16) 15 111110010110000010101111100100 1.0460E*09 0.7694
1 15> 111110011010000010011111110001 1.0470E+09 0.7772 15) (20, 1) 30 111111001000000010010011000101 1.0591E*09 0.8715
j 16) 111111001010000010101111100100 1.0596E+09 0.8758 16) (20, 1) 30 111100011010010010101111110001 1.0135E*09 0.5615
i 17) 111110011110000010110011101111 1.0481E+09 0.7850 17) ( 9,10) 30 111111101000000010011011100111 1 .0675E*09 0.9430
1 18) 111111011000001110101111100100 1.0633E+09 0.9070 18) ( 9,10) 30 111111001100000010011111110111 1.0601E*09 0.8801
19) 111111001000000010011111110000 1.0591E+09 0.8715 19) ( 3, 8) 14 111111011000000010011111110100 1.0633E+09 0.9066
J 20) 111111001000000010010011000101 1.0591E+09 0.8715 20) ( 3, 8) 14 111111001000000010011011000100 1.0591E+09 0.8715
J 21) 111110011110000010101111110000 1.0481E+09 0.7850 21) { 1,26) 8 111110010110000000001110110001 1.0460E*09 0.7694
22) 101110011110010010101111100100 7.7969E+08 0.0408 22) ( 1,26) 8 111100011010010010101111111111 1.0135E*09 0.5615
23) 111111101000000010011111110001 1.0675E*09 0.9430 23) (18,20) 30 111011011000001110101111110100 9.9621E*08 0.4726
24) 111111001000000010011111110111 1.0591E+09 0.8715 24) (18,20) 30 111111001000000010010010000101 1.0591E+09 0.8715
25, 111111111000000010011111100100 1.0717E+09 0.9807 25) (14,30) 3 111111001001100000101111100000 1.0595E+09 0.8747
26) 111110010110000000101111110111 1.0460E+09 0.7694 26) (14,30) 3 111111001001110000101111010100 1.0595E*09 0.8752
27) 111111001001110001101111110000 1.0595E*09 0.8752 27) (23, 3) 18 111111001001000010011111110001 1.0593E+09 0.8736
28) 011110011000010010000011000101 5.0968E-i-08 0.0006 28) (23, 3) 18 111111101100000010011111110100 1.0685E*09 0.9523
: 111111001100000010011011110001 1.0601E+09 0.8801 29) (24, 2) 30 111111001000010010001111110111 1.0591E+09 0.8720
30) i 111111001001100000101111100100 1.0595E*09 0.8747 30) (24, 2) 30 111111001000100010011111110011 1.0592E*09 0.8726
| Notę: Generation 7 & Accimulated Statistics: max=0.9807, min=0.0405, avg=0.8100, sum=24.2997, nmutation�201, ncross* 71
Wyd. 3.14. SGA, raport o pokoleniach 6 i 7
90
3. Implementacja komputerowa algorytmu genetycznego
zbieżność płacą za to utratą globalności i giętkości. Wiele z tych metod ma ograniczone pole zastosowań. W konsekwencji algorytmy genetyczne nadają się do stosowania tam, gdzie metody zbieżne nie odważą się nawet stąpnąć. Co więcej, w przypadku rozwiązywania zadań, dla których są znane zbieżne, ale lokalne algorytmy, można w naturalny sposób myśleć o hybrydyzacji. Należy rozpocząć poszukiwania przy użyciu algorytmu genetycznego, aby odszukać interesujące wzniesienia terenu. Gdy algorytm genetyczny wywęszyjuż, które rejony zapowiadają się najlepiej, należy wyciągnąć lokalny zbieżny algorytm i wykonać wspinaczkę na miejscowe wierzchołki. W ten sposób można połączyć globalność i równoległość algorytmu genetycznego ze zbieżnością techniki lokalnej. Za chwilę zajmiemy się jednym z praktycznych sposobów redukcji przedwczesnej zbieżności za pomocą skalowania funkcji przystosowania. Najpierw jednak musimy ob-myśleć technikę przekształcania dowolnej funkcji celu do postaci wymaganej dla funkcji przystosowania. �.� > � �
3.6. Przekształcenie funkcji celu w funkcję przystosowania _
Wiele zadań wyraża się w bardziej naturalny sposób w kategoriach minimalizacji pewnej funkcji kosztu g(x) niż maksymalizacji pewnej funkcji użyteczności lub zysku u(x). Nawet w przypadku, gdy mamy w naturalny sposób do czynienia z zadaniem maksymalizacji, nie gwarantuje to, że funkcja użyteczności będzie przyjmowała wartości nieujemne dla wszystkich jc, co musi zachodzić dla funkcji przystosowania ^ak pamiętamy, funkcja przystosowania jest nieujemnym kryterium jakości). Wskutek tego często zachodzi potrzeba przekształcenia oryginalnej funkcji celu w funkcję przystosowania, czego można dokonać w jednym lub kilku krokach.
Dualność zadań minimalizacji kosztu i maksymalizacji zysku jest dobrze znana. W badaniach operacyjnych, chcąc przejść od zadania minimalizacji do zadania maksymalizacji, mnożymy po prostu funkcję kosztu przez minus jeden. W naszym przypadku sama taka operacjajest niewystarczająca, gdyż otrzymany miernik nie musi być zawsze nieujemny. Na użytek algorytmów genetycznych stosuje się więc powszechnie następujące przekształcenie kosztu w przystosowanie:
/(*) =
' Cmta - g(x), jeżeli g(x) < Cmax
0
w przeciwnym przypadku
Współczynnik Cmm można wybierać na różne sposoby. Może on być wczytany z wejścia, może być równy największej napotkanej do tej pory wartości funkcji g, albo największej wartości g w populacji bieżącej lub w k ostatnich populacjach. Być może właściwiej byłoby, żeby CmaK wahał się w zależności od wariancji populacji. Tę ostatnią możliwość rozważymy w rozdziale czwartym.
Kiedy w zadaniu występuje w naturalny sposób funkcja zysku lub użyteczności, nie mamy kłopotów z "kierunkiem" funkcji: maksymalizacja zysku lub użyteczności
3.7. Skalowanie przystosowania
91
prowadzi do pożądanego celu. Możemy mieć jednak problemy z ujemnymi wartościami funkcji użyteczności u(x). Trudności te łatwo można pokonać za pomocą przekształcenia
f(x) = \
0
+ Caln, jeżeli u(x) + Cm(, > 0
w przeciwnym przypadku
Współczynnik Cmin może być wczytany z wejścia, może być równy modułowi najmniejszej wartości u w populacji bieżącej lub ostatnich k populacjach, wreszcie może być funkcją wariancji populacji. Dalszą dyskusję tych możliwości odłożymy do następnego rozdziału.
Wszystkie te manipulacje z funkcjami celu mogą zrodzić podejrzenia co do istoty związku między funkcjami celu a funkcjami przystosowania. W przyrodzie przystosowanie (liczba potomstwa dorastającego do wieku reprodukcyjnego) jest tautologią. Wiele potomstwa przeżywa, gdyż rodzice są dobrze przystosowani, a rodzice są dobrze przystosowani, gdyż wiele potomstwa przeżywa0. W naturalnych populacjach przeżycie jest ostatecznym i jedynym kryterium oceny wszelkich nowinek. W przeciwieństwie do tego, podczas działania algorytmu genetycznego mamy możliwość, a może i obowiązek regulowania poziomu współzawodnictwa między członkami populacji w celu osiągnięcia pożądanego zachowania w fazach przejściowej i końcowej procesu. Dokładnie o to właśnie chodzi, gdy dokonujemy skalowania przystosowania.
3.7. Skalowanie przystosowania
Regulacja liczby kopii ma szczególne znaczenie dla algorytmów genetycznych pracujących na małych populacjach. Często się zdarza, że na początku przebiegu w populacji występuje kilku ponadprzeciętnych osobników, podczas gdy reszta zalicza się do pośledniej kategorii. Gdyby zastosować normalną regułę selekcji &select,=f|/^f), te ponad-przeciętne osobniki uzyskałyby znaczący udział w skończonej populacji w ciągu jednego pokolenia, a to jest zjawisko niepożądane, główna przyczyna przedwczesnej zbieżności. W późniejszej fazie mamy do czynienia z zupełnie odmiennym problemem. Może się zdarzyć, że pod koniec przebiegu populacja zachowała znaczną różnorodność, ale średni wskaźnik przystosowania niewiele odbiega od maksymalnego. Jeśli na to nie zareagujemy, to osobniki przeciętne i osobniki najlepsze będą otrzymywać prawie tę samą liczbę potomstwa w następnych pokoleniach, co sprowadzi zasadę przeżycia najlepiej przystosowanych siłę sprawczą ulepszeń - do poziomu błądzenia przypadkowego wśród przeciętniaków. Skalowanie przystosowania może być pomocne w obu sytuacjach - na początku i w fazie dojrzałości.
" Pojęcie "najlepiej przystosowanego" (pochodzące od Spencera) natrafiało wielokrotnie na trudności interpretacyjne. Współczesne naświetlenie tego zagadnienia można znaleźć np. w książce E. Mayra Populacje, gatunki i ewolucja, WP, Warszawa l974, str. 150 tyrzyp. tlum.).
92
3. Implementacja komputerowa algorytmu genetycznego
"f Jedną z użytecznych metod skalowaniajest skalowanie liniowe. Niech/oznacza przystosowanie pierwotne, a/' - przystosowanie po skalowaniu. W przypadku skalowania liniowego, związek miedzy / a /' musi mieć postać
f' = af+b
Współczynniki a i b można wybrać na wiele sposobów; w każdym razie chcemy jednak, aby średnie przystosowanie po skalowaniu/"V4, było równe pierwotnemu, tj.fmg, gdyż zapewnia to, że osobnik o przeciętnym przystosowaniu będzie miał średnio jednego potomka w następnym pokoleniu. Średnią liczbę potomków osobnika o maksymalnym (pierwotnym) przystosowaniu możemy kontrolować za pomocą drugiego warunku
f' - C � f � �� ��
J mas. ^~mult Javg
gdzie Cmult (współczynnik zwielokrotnienia) jest właśnie żądaną liczbą kopii. W typowych, małych populacjach (od 50 do 100 osobników) dobre wyniki otrzymywano, wybierając wartość Cmull między 1,2 a 2,0.
, ...,,,,- 'min 'avg 'max
Przystosowaniepierwotne ,
Rys. 3.5. Skalowanie liniowe w typowym przypadku
Taki dobór Cmull powoduje znaczne rozciągnięcie funkcji przystosowania pod koniec przebiegu. Może to spowodować pewne kłopoty z zastosowaniem reguły skalowania liniowego, zilustrowanej na rys. 3.5. Początkowo nie ma problemów z przekształceniem, gdyż kilku ponadprzeciętnych osobników spada w dół, a gorsi osobnicy wędrują nieco w górę. Z trudniejszym przypadkiem mamy do czynienia na rys. 3.6. Takie sytuacje są typowe w dojrzałych populacjach, w których nieliczni "degeneraci" (ciągi o kiepskim przystosowaniu) wypadają znacznie poniżej średniej, niewiele odbiegającej od maksimum. Gdy zastosujemy w takim przypadku regułę skalowania liniowego, rozciągając w wymaganym stopniu funkcję przystosowania, małe wartości funkcji stają się ujemne.
3.7. Skalowanie przystosowania
93
2f'
c 0!
avg
'avg
� f'r
avg 'max Przystosowanie pierwotne
Ujemne wartości przystosowania naruszająwymóg nieujemności
Rys. 3.6. Efekty skalowania liniowego w dojrzałej populacji. Punkty o małym przystosowaniu mogą przejść w obszar wartości ujemnych
Możemy sobie z tym różnie radzić; tutaj przyjmujemy, żejeśli skalowanie ze współczynnikiem zwielokrotnienia Cmull okazuje się niewykonalne, zachowujemy nadal równość obu średnich, ajako drugi warunek wybieramy równość/^," = 0 (minimalne przystosowanie pierwotne przy skalowaniu staje się równe zeru).
Programbędący implementacją elementarnego algorytmu genetycznego można bez trudu uzupełnić o skalowanie, korzystając z trzech prostych podprogramów prescale, scale oraz scalepop (wyd. 3.15). Procedurapresca/e oblicza na podstawie wskaźników przystosowania średniego, maksymalnego i minimalnego (nazwanych odpowiednio uavg, umax i umin) współczynniki przekształcenia liniowego, a i b, zgodnie z opisanymi wyżej regułami. Najpierw sprawdza się, czy wykonalne jest skalowanie z pożądanym współczynnikiem zwielokrotnienia Cmult (występującym w programie pod nazwą finu-ltiple)\ jeśli tak, to obliczenia idą odpowiednim torem. W przeciwnym przypadku skalowanie wykonuje się, zachowując średnią i "rozciągając" przystosowanie w obie strony - tak by minimum przeszło na zero. Po wywołaniu przygotowawczej procedury prescale przychodzi kolej na procedurę scalepop, która, za pośrednictwem prostej funkcji paskalowej scale, dokonuje przeskalowania wskaźników przystosowania wszystkich osobników w populacji. Zakładamy tu, że pierwotne wskaźniki przystosowania są zapisane w polu objective rekordów typu individual, reprezentujących osobników (zmienne pop^].objective). Zmodyfikowane wskaźniki przystosowania zostają umieszczone w polu fitness tych samych rekordów, a przy okazji wyznacza się nową wartość łącznej sumy wskaźników sumfitness. Instalację i przetestowanie procedur skalowania pozostawiamy jako ćwiczenie.
W taki oto sposób proste skalowanie pomaga zapobiegać przedwczesnej dominacji nadmiernie wybijających się osobników, następnie zaś pobudza do zdrowej konkurencji tych o wyrównanych szansach. Nie zamykamy jeszcze naszej dyskusji na temat możliwych przekształceń funkcji celu; wrócimy do niej w rozdziale czwartym, podając kilka nowych przykładów. A teraz zajmiemy się kilku sposobami kodowania, które nadają się do zastosowania w algorytmach genetycznych - innymi niż używany do tej pory prosty kod.
94
3. Implementacja komputerowa algorytmu genetycznego
{ scale.sga: contains prescale, scale, scalepop for scaling fitnesses )
procedurę prescale(umax, uavg, umin:real; var a, b:real); ( Calculate scaling coefficients for linear scaling ) const fmultiple - 2.0; ( Fitness multiple is 2 ) var delta:real; ( Divisor ) begin
if umin > (fmultiple*uavg - umax) / (fmultiple - 1.0) ( Non-negative test ) then begin ( Normal Scaling ) delta : umax - uavg;
a :- (fmultiple - 1.0) * uavg / delta; ;
b : uavg * (umax - fmultiple*uavg) / delta; end else begin ( Scale as much as possible ) delta :- uavg - umin;
a
b
end
- uavg / delta;
-umin * uavg / delta;
end;
function scale(u, a, b:real);real; ( Scale an objective function value ) begin scale :- a * u + b end;
procedurę scalepop(popsize:integer; var max, avg, min, sumfitness:real;
var pop:population); ( Scale entire population ) .... ,,. var j :integer;
a, b:real; ( slope & intercept for linear equation ) begin
prescale(max, avg, min, a, b); ( Get slope and intercept for function ) sumfitness :- 0.0;
for j :- 1 to popsize do with pop[j] do begin fitness :- scale(objective, a, b); sumfitness :- sumfitness + fitness; end; end;
Wyd. 3.15. Podprogramy skalowania: procedura prescale, funkcja scale i procedura scalepop
3.8. Metody kodowania
Jak dotąd mieliśmy do czynienia z bardzo ograniczonym repertuarem sposobów odwzorowania ciągów kodowych na parametry zadań optymalizacyjnych. W zadaniu o przełącznikach użyliśmy prostej metody kodowania binarnego, w której zera lub jedynki na poszczególnych pozycjach ciągu kodowego odpowiadały położeniom odpowiednich przełączników (wyłączony/włączony). W drugim przypadku interpretowaliśmy ciąg zer i jedynek jako nieujemną liczbę całkowitą, przy czym ciągowi A = a,a,_\ ... a2at odpowiadała wartość parametru x = ^a^2'^ (zapis odwrócony!). Chociaż kody te dają pewną swobodę, nie są wystarczająco giętkie, aby poradzić sobie z całym spektrum problemów, z którymi spotykamy się w nauce, technice i interesach. W tym punkcie omówimy dwie podstawowe zasady kodowania na użytek algorytmów genetycznych, co powinno w przyszłości ułatwić projektowanie kodów dla różnych typów zadań. Nieco później przedstawimy standardową metodę kodowania binarnego dla zadań wielopara-
3.8. Metody kodowania
95
metrycznych, która dowiodła - i powinna nadal dowodzić - swej użyteczności na przykładzie rozmaitych zadań.
W pewnym sensie wybór kodu dla zadania rozwiązywanego przy użyciu algorytmu genetycznego nie stanowi żadnego problemu, gdyż programista jest ograniczony głównie swoją własną wyobraźnią. Z poprzedniego rozdziału wiemy, że algorytmy genetyczne potrafią wyzyskać podobieństwa ciągów kodowych - niezależnie od metody kodowania -jeśli tylko ze schematów-cegiełek (tj. wąskich schematów o wysokim przystosowaniu) można budować niemal optymalne rozwiązania. Patrząc z innego punktu widzenia, ta swoboda wyboru stanowi wątpliwe błogosławieństwo dla niedoświadczonego użytkownika; widok całych szeregów możliwych sposobów kodowania może zarówno dodawać otuchy, jak i oszałamiać. Mając taką wolność wyboru - jak nowy użytkownik będzie w stanie znaleźć dobry kod? Na szczęście, algorytmy genetyczne są pod tym względem dość wyrozumiałe, ponieważ są odporne; dlatego na ogół nie ma sensu dręczyć się problemem wyboru metody kodowania. Na dodatek, proponujemy kierować się dwiema podstawowymi zasadami konstrukcji kodów dla algorytmów genetycznych. Są to: zasada znaczących cegielek \ zasada minimalnego c.lfabetu.
Zasada znaczących cegiełek stwierdza, co następuje: ...... .
Kod należy dobierać w taki sposób, żeby schematy niskiego rzędu i o małej rozpiętości wyrażaty własności zadania oraz pozostawaty względnie niezależne od schematówoinnychpozycjachustalonych.
Choć dzięki analizie Walsha, o której wspomnieliśmy w rozdziale drugim, można sprawdzić, czy dany kod czyni zadość powyższej zasadzie, podejście tojest niezbyt praktyczne, tak więc projektowanie kodu zapewniającego znaczące cegiełki jest czymś w rodzaju sztuki. Niemniej jednak projektując kod powinniśmy zwrócić uwagę na dystans między powiązanymi logicznie pozycjami. W rozdziale piątym omawiamy sposoby rekonfigura-cji ciągów kodowych, jak również kilka operacji, dla których poszukiwanie dobrych rozwiązań sprowadza się do poszukiwania dobrych kodów.
Druga zasada projektowania kodów, zasada minimalnego alfabetu, mówi po prostu: ; :
Należy wybrać najmniejszy alfabet, w którym dane zadanie wyraża się w sposób naturalny.
Do tej pory trzymaliśmy się niemal obsesyjnie idei kodu dwójkowego. Czy był to tylko przypadek, czy też w naszym kodowym szaleństwie była metoda? Że była to jednak metoda, najlepiej możemy się przekonać na wyeksploatowanym już, ale pouczającym przykładzie kodu pięciopozycyjnego, który wprowadziliśmy w rozdziale pierwszym. W tablicy 3.1 widzimy cztery znajome dwójkowe ciągi kodowe o tych samych wskaźnikach przystosowania, co poprzednio. Jak pamiętamy, jednym z naszych pierwotnych motywów wprowadzenia pojęcia schematu była zrozumiała chęć powiązania wysokiego przystosowania ze strukturą schematów występujących w populacji. We wspomnianej tablicy zamieściliśmy też pewien kod niedwójkowy. W istocie, chodziło tu o przypadek
96
3. Implementacja komputerowa algorytmu genetycznego
krańcowy. Rozważmy wzajemnie jednoznaczne przyporządkowanie pięcioelementowych ciągów zerojedynkowych i symboli alfabetu złożonego z 26 liter A-Z oraz sześciu cyfr 1-6, jak pokazano w tabl. 3.2.
Tablica 3.1. Porównanie ciągów kodowych dla kodu dwójkowego i niedwójkowego
Ciąg dwójkowy
Wartość X Ciąg niedwójkowy Przystosowanie
01101 13 N 169
11000 24 Y 576
01000 8 I 64
10011 19 T 361
Przeglądając listę ciągów (tabl. 3.1) możemy zauważyć, że w przypadku zerojedyn-kowym wykrycie podobieństw strukturalnych stało się możliwe dzięki małemu alfabetowi. Natomiast w drugim przypadku mamy do czynienia z czterema pojedynczymi literami (ciągami jednoelementowymi) i odpowiadającymi im wskaźnikami przystosowania; nie widać tu żadnych podobieństw, które moglibyśmy wykorzystać. Prawda, przykład jest krańcowy, ale ta sama zasada obowiązuje w mniej jaskrawych przypadkach.
Tablica 3.2. Tabela odpowiedniości dla kodu dwójkowego i niedwójkowego
Tabela odpowiedniości kodów
dwójkowy
niedwójkowy
00000 00001
A B
11001 11010 11011
z 1
2
11111
Patrząc na to od strony matematycznej, spróbujmy porównać liczbę schematów występujących w przypadku kodów dwójkowych i niedwójkowych. Oczywiście, w obu przypadkach powinniśmy mieć tę samą liczbę ciągów kodowych; jednakże różnica wielkości alfabetów pociąga za sobą różną długość ciągów kodowych. Aby liczba punktów w obu przestrzeniach byłajednakowa, musi zachodzić warunek 2' = fc', gdzie /jest długością ciągu zerojedynkowego, a /' - długością słowa w alfabecie fc-elementowym. Liczbę schematów dla obu kodów można wyznaczyć na podstawie długości ciągu: dla alfabetu dwuelementowego wynosi ona 3', dla &-elementowego zaś (&+!)'. Nietrudno udowodnić, że alfabet dwuelementowy charakteryzuje się największą ze wszystkich licz-
3.9. Standardowa metoda kodowania dla zadań wieloparametrycznych
97
bą schematów przypadających na bit informacji. Ponieważ owe schematy leżą u sedna procesu poszukiwań, więc projektując kod powinniśmy starać się zmaksymalizować ich liczbę ~ w ten sposób algorytm genetyczny będzie miał więcej materiału do eksploatacji.
3.9. Standardowa metoda kodowania dla zadań ____________________ wieloparametrycznych
Sformułowane w p. 3.8 dwie zasady dostarczają ogólnych wskazówek na temat projektowania efektywnych kodów dla algorytmu genetycznego, nie określają jednak praktycznych sposobów doboru kodu dla konkretnego zadania. Jedną ze skutecznych metod kodowania zadań optymalizacyjnych z wieloma parametrami jest użycie blokowego zapisu pozycyjnego ze standaryzacją parametrów.
Mieliśmy już do czynienia z kodowaniem nieujemnych liczb całkowitych z przedziału [0, 2'- 1], cojednak zrobić w innych przypadkach? Jednym ze sposobów pozbycia się tego ograniczenia jest liniowe odwzorowanie zdekodowanej liczby z przedziału [0, 2'~1] w zadany przedział [Umin, Umax\. Dzięki temu jesteśmy w stanie starannie kontrolować zakres i precyzję zmiennych decyzyjnych. Dokładność opisanego kodu jest równa
TC =
JI U
max min
2'-l
Aby skonstruować kod w przypadku wieloparametrycznym, wystarczy po prostu skon-katenować bloki kodowe reprezentujące poszczególne parametry. Każdy blok może mieć, oczywiście, inną długość i reprezentować inny przedział [Umin, Umm] (rys. 3.7).
Zestaw podprogramów implementujących blokowy zapis pozycyjny ze standaryzacją parametrów został podany w wyd. 3.16. Korzystają one z opisanej wcześniej funkcji decode w celu zdekodowania do postaci standaryzowanej pojedynczego bloku kodowego reprezentującego jeden parametr. Procedura extract-parm wydobywa pojedynczy blok
Jeden parametr, przedział L/, (/, = 4)
_____________________________ .Nf*ii >� '��
0 0 0 0 -------* Unin . v, . .-� ��: -^
1 1 1 1 -------* Unax
Pozostałe ciągi przechodzą liniowo na wartości pośrednie
Kod blokowy (10 parametrów)
0 0 0 1 ] 0 1 0 1 ]
I
] 1 1 0 0 j 1 1 1 1 J
I
U10
I
I I I I " I
Rys. 3.7. Blokowy zapis pozycyjny ze standaryzacją dla zadań wieloparametrycznych
98
3. Implementacja komputerowa algorytmu genetycznego
z pełnego ciągu kodowego, funkcja paskalowa map-parm dokonuje przekształcenia wartości standaryzowanej na liczbę z przedziału [minparm, maxparm], a procedura deco-de~parms koordynuje proces dekodowania wszystkich nparms parametrów. Instalację i testowanie tych podprogramów pozostawiamy jako ćwiczenie.
type parmparra - record ( parameters of the pararaeter ) lparm:integer; ( length of the parameter ) parameter, maxparm, minparm:real; ( parameter & range ) end; parraspecs array[l..maxparms] of parmparm;
var parras:parmspecs; � - . - .., ; � -.: . ; -� ., : �. y. , � \ .;� . ����, < �.-, -;
; procedurę extract_parm(var chromfrom, chromto:chroraosome;
var jposition, lchrom, lparm:integer); ( Extract a substring from a full strlng )
var j, jtarget:integer; �; j ,� '
begin
j :- l; � ' " '
jtarget :- jposition + lparm - 1; .., if jtarget > lchrom then jtarget :- lchrora; ( Clarap if excessive )
while (jposition <- jtarget) do begin
'? chromto[j] :- chromfrom[jposition]; � �>'.�"�
jposition :- jposition + 1; j :- j + 1; end; end;
functionmap_parm(x, maxparm, minparm, fullscale:real):real; { Map an unsigned binary integer to range [minparm,maxparmJ ) begih map_parm :- minparm + (maxparm - minparm)/fullscale*x end;
procedurę decode_parms(var nparms, lchrom:integer; var chrom:chromosome; var parms:parmspecs); var j, jposition:integer; ,;
chromtemp:chromosome; ( Temporary string buffer ) begin
j : 1; ( Parameter counter ) ' � ,*
i , jposition :- 1; ( String position counter } , . lt
repeat ' ' " ' '
with parms[j] do if lparm>0 then begin
extract_parm(chrom, chromtemp, jposition, lchrom, lparm); parameter :-map_parm( decode(chromtemp, lparm),
maxparm, minparm, power(2.0, lparm)-1.0 ); end else parameter :- 0.0;
J :-j + 1; .,,-.�,*.. .-,,
until j > nparms; end;
Wyd. 3.16. Podprogramy dekodujące do zastosowania w SGA: procedura extracLparm, funkcja map_pam? i procedura decode^>arms
3.10. Dyskretyzacja
99
3.10. Dyskretyzacja
Dyskretyzacja parametrycznego zadania optymalizacji z parametrami rzeczywistymi nie jest jedynym rodzajem dyskretyzacji, która może być konieczna w celu zastosowania algorytmu genetycznego. Wiele zadań optymalizacyjnych, a właściwie zadań z dziedziny optymalnego sterowania, charakteryzuje się nie pojedynczym parametrem decyzyjnym, ale funkcją decyzyjną (zwaną sterowaniem), której wartość musi być określona w każdym punkcie pewnego kontinuum. Zastosowanie algorytmów genetycznych do rozwiązywania tego typu problemów jest możliwe po uprzednim sprowadzeniu ich do postaci parametrycznej (skończenie wymiarowej).
Taką formę dyskretyzacji można prosto zilustrować następującym przykładem: Przypuśćmy, że naszym celem jest minimalizacja czasu przejazdu rowerem danego odcinka drogi. Załóżmy dalej, że możemy przykładać siłę/jako funkcję czasu/(f), przy ograniczeniu I f(f} \^fmax. Jest to ciągłe zadanie optymalnego sterowania, w którym chodzi o znalezienie wielkości przykładanej siły, wyrażonej jako ciągła funkcja czasu (por. rys. 3.8). Ponieważ stosowanie algorytmu genetycznego jest możliwe tylko dla zadań dających się opisać za pomocą struktur skończonych, musimy najpierw sprowadzić problem ciągły do skończonej liczby parametrów, a następnie zredukować te parametry do postaci skończonych ciągów kodowych przy użyciu pewnego procesu kodowania.
Stta F
Rozwiązanie ciągłe
to
t,
t.
t,
Czas t ...... ,...
Rys. 3.8. Dyskretyzacja ciągłej funkcji decyzyjnej
Dyskretyzacja kontinuum, o której tu mówimy, kojarzy się najczęściej z takimi tematami, jak sterowanie dyskretne, interpolacja i metoda elementu skończonego. W zadaniu sterowania rowerem pierwszy etap dyskretyzacji mógłby polegać na zastąpieniu funkcji ciągłej ciągiemjej wartości/j, wziętych w równych odstępach czasu. Następnie możemy wybrać pewną postać funkcyjną, np. schodkową, przedziałami liniową, przedziałami kwadratową albo splajn trzeciego stopnia i przeprowadzić interpolację na podstawie wartosci_^. Na rysunku 3.8 pokazano przybliżone rozwiązanie zadania sterowania rowerem będące interpolacją liniową rozwiązania dokładnego.
100
3.11. Problem więzów
3. Implementacja komputerowa algorytmu genetycznego
Do chwili obecnej ograniczyliśmy się jedynie do rozpatrywania bezwarunkowych zadań optymalizacyjnych. W praktycznych problemach często występujejedno lub więcej ograniczeń, które muszą być spełnione. W tym punkcie zajmiemy się sprawą uwzględnienia takich warunków (zwanych więzami) w kontekście zastosowania algorytmu genetycznego.
Więzy mogą przybierać postać równości lub nierówności. Ponieważ więzy o charakterze równości można włączyć do modelu systemu - czarnej skrzynki - pozostaje nam rozpatrzyć tylko więzy wyrażone w postaci nierówności. Na pozór wydaje się, że takie ograniczenia nie powinny przedstawiać szczególnego problemu. Algorytm genetyczny generuje ciąg parametrów, które są poddawane ocenie, biorąc pod uwagę model systemu, funkcję celu i więzy. Należy więc każdorazowo obliczyć wartość funkcji celu i sprawdzić, czy nie zostały naruszone któreś z więzów. Jeśli nie, to przypisujemy danemu zbiorowi parametrów wskaźnik przystosowania odpowiadający wyznaczonej wartości funkcji celu. Jeśli natomiast więzy zostały naruszone, to dane rozwiązanie jest niedopuszczalne, nie można mu więc przypisać przystosowania. Byłaby to dobra reguła postępowania, gdyby nie to, że sporo zadań spotykanych w praktyce to zadania z silnymi więzami; znalezienie rozwiązania dopuszczalnego może być dla nich niemal tak samo trudne,jak znalezienie rozwiązania optymalnego. W rezultacie często zmuszenijesteśmy wydobywać jakąś informację z rozwiązań niedopuszczalnych, być może obniżając ich wskaźnik przystosowania w stopniu odpowiadającym naruszeniu więzów. Tak postępuje sie,stosu)acmetodekarfy>enaltymethod]. .�'
W metodzie kar zadanie optymalizacyjne z więzami przekształcamy do postaci bezwarunkowej, obciążając każde naruszenie więzów kosztem lub karą. Koszt ten zostaje włączony do funkcji celu. Rozważmy dla przykładu następujące warunkowe zadanie minimalizacji: ,
zminimalizować g(x)
przy warunku h|(x) > 0, /=1, 2, ..., n ',
gdzie je jest wektorem m-wymiarowym. ;
<.-:.A
Po przekształceniu do postaci bezwarunkowej:
n
zminimalizować g(x) + rX<&[fy(;c)]
/=i gdzie O - funkcja kary
r - współczynnik kary.
Funkcję kary O można wybrać na wiele sposobów. W tej książce najczęściej używamy funkcji kwadratowej ^>[h/(x)] = h^(x), dla wszystkich naruszonych więzów /. Pod pewnymi warunkami ciąg rozwiązań zadania bezwarunkowego jest zbieżny do rozwiązania zadania warunkowego, gdy współczynnik kary r dąży do nieskończoności. W praktyce algorytmów genetycznych wartości r ustala się odrębnie dla każdego typu więzów, tak aby umiarkowane przekroczenia więzów były obciążone karą stanowiącą ustalony, znaczący ułamek kosztu nominalnego. " ..> >
3.12. Podsumowanie
101
3.12. Podsumowanie
W rozdziale tym uchyliśmy nieco zasłonę tajemnicy algorytmów genetycznych, dzięki uważnemu zbadaniu struktur danych, procedur i rozmaitych szczegółów niezbędnych do implementacji prostego, ale dającego się praktycznie stosować algorytmu. W szczególności zaprogramowaliśmy w Pascalu szkielet algorytmu genetycznego w postaci programu SGA, który nadaje się do wykonania przy użyciu powszechnie dostępnych mikrokomputerów.
Jak można się było spodziewać, podstawową strukturą danych dla elementarnego algorytmu genetycznego jest populacja ciągów kodowych. W programie SGA używa się dwóch identycznych struktur, aby w maksymalnym stopniu ułatwić tworzenie i wymianę populacji. Sama populacja jest zrealizowana w postaci tablicy osobników, złożonych z ciągu kodowego, zdekodowanego parametru, wskaźnika przystosowania i pewnych ważnych informacji pomocniczych.
Zasadnicza część pracy programu przypada na trzy podprogramy: select, crossover i mutation. Pierwszy z nich dokonuje wyboru losowego z powtórzeniami zgodnie z zasadą zwaną regułą ruletki. Dwa następne realizują, zgodnie ze swymi nazwami, operacje opisane w rozdziale pierwszym. Wszystkie działania są koordynowane przez procedurę o nazwie generation, która w każdym kolejnym pokoleniu generuje nową populację.
Po skonstruowaniu SGA przetestowaliśmy program w krótkim przebiegu, dla nieskomplikowanej funkcji/(x) = cx10. Chociaż nie można wyciągać ogólnych wniosków na podstawie jednej realizacji procesu stochastycznego, algorytm znalazł rozwiązanie bliskie optymalnemu w krótkim czasie i po przeszukaniu drobnej cząstki przestrzeni poszukiwań.
Omówiliśmy również rozmaite szczegóły implementacyjne, w tym przekształcenie funkcji celu do postaci funkcji przystosowania, które trzeba niekiedy wykonać dla tradycyjnie sformułowanych zadań optymalizacyjnych (np. w przypadku zadań minimalizacji), aby nadawały się do rozwiązania za pomocą algorytmu genetycznego. Ponadto zaproponowaliśmy skalowanie wskaźników przystosowania w celu dokładniejszej kon-trolitempareprodukcjinajlepszychciągówkodowych.
Przedyskutowaliśmy też pewne podstawowe kwestie dotyczące problemu kodowania. Zostały sformułowane zasady minimalnego alfabetu \ znaczących cegielek, określające kryteria projektowania efektywnych kodów. Zaprezentowaliśmy zaprogramowaną w Pascalu metodę kodowania, przydatną w rozmaitych zadaniach optymalizacyjnych. W metodzie tej używa się blokowego zapisu pozycyjnego ze standaryzowaniem parametrów.
Kodowanie pociąga za sobą dyskretyzację parametrów, ale często niejest tojedyny rodzaj dyskretyzacji, z jakim możemy się spotkać. W licznych zadaniach optymalizacji mamy do czynienia z ciągłą funkcją decyzyjną. Takie zadania optymalizacji - ściślej, zadania optymalnego sterowania - trzeba najpierw sprowadzić do postaci skończenie wymiarowej. Otrzymane w ten sposób parametry mogą być następnie zakodowane metodą opisaną w tym rozdziale. Dyskretyzację tego typu wykonuje się wybierając odpowiednie funkcje interpolacyjne (zazwyczaj wystarcza interpolacja liniowa); parametry związane
102
3. Implementacja komputerowa algorytmu genetycznego
z wybraną postacią funkcji interpolacyjnej po zakodowaniu przy użyciu kodu blokowego tworzą następnie skończony dwójkowy ciąg kodowy.
Na zakończenie uznaliśmy potrzebę zastosowania specjalnych metod w celu uwzględnienia więzów wyrażonych w postaci nierówności. Naturalną dziedziną zastosowania algorytmów genetycznych są problemy bezwarunkowe. Wszak przyroda bezustannie eksperymentuje, a jedynym rodzajem więzów, jakie napotyka w swym środowisku jest kryterium przetrwania jej wytworów. Także w przypadku algorytmów genetycznych funkcja przystosowania musi odzwierciedlać zarówno jakość, jak i dopuszczalność rozwiązania. W licznych przypadkach skuteczna okazała się tu metoda kar. Przy użyciu tej metody każde naruszenie więzów jest karane przez obciążenie pierwotnej funkcji celu kosztem zależnym od miary przekroczenia.
Analiza konkretnej implementacji i dyskusja niektórych zagadnień implementacyjnych zbliżyły nas do praktycznego zastosowania algorytmów genetycznych w nauce, technice i interesach. W następnym rozdziale przyjrzymy się wcześniejszym i obecnym przykładom takich zastosowań.
3.13. Zadania
3.1. Rozważmy funkcję przystosowania/(*)=*" w przedziale [0, 1]. Oblicz wartość oczekiwaną średniego przystosowania dla populacji losowo wygenerowanych punktów. Oblicz prawdopodobieństwo wygenerowania punktu x, dla którego f(x) >f. Porównaj wartości liczbowe średnich z populacji dla n = 2 i n= 10. Porównaj prawdopodobieństwa wygenerowania punktu o wskaźniku przystosowania/>0,9 w obu tychprzypadkach. ;;,�;;;' ^ ,-
3.2. Przestrzeń rozwiązań zawiera 2097152 punkty. Chcemy porównać dwie wersje algorytmu genetycznego, jedną z kodem dwójkowym i drugą z kodem ósemkowym. Oblicz i porównaj następujące wielkości dla obu wersji algorytmu:
� : � . l st ; � ' *
a) całkowitaliczbaschematów; ' '-* �
b) całkowita liczba ciągów kodowych;
c) liczba schematów przypadających na jeden ciąg kodowy;
d) ograniczenie dolne i górne liczby schematów reprezentowanych w populacji liczącej 50 elementów.
; �. " . iWrffiłt f'>iffi'f ':f'-. � ' ,� � '
3.3. Chcemy zminimalizować funkcję trzech zmiennych,/(*,y,z). Zmienna x może przybierać wartości z przedziału [-20,0, 125,0], zmienna y - z przedziału [0, l,2x 106], a zmienna z - z przedziału [-1,0, l,0]. Wymagana dokładność dla x, y \ z wynosi odpowiednio 0,5, 104 i 0,001. Dobierz odpowiedni dla tego zadania blokowy zapis pozycyjny ze standaryzacją parametrów. Jakajest minimalna liczba bitów niezbędna do osiągnięcia żądanej dokładności? Podaj ciągi kodowe reprezentujące następujące punkty: (-20, 0, -1), (125,0, l,2xl06, 1,0), (50, 100000, 0,597).
3.14. Ćwiczenia komputerowe
103
3.14. Ćwiczenia komputerowe
A. Włącz do programu SGA podprogamy skalujące z wyd. 3.15 i powtórz eksperyment z funkcja/(*) = c*'". Porównaj i omów wyniki otrzymane w obu wariantach (bez skalowania i ze skalowaniem).
B. Przetestuj podprogramy realizujące blokowy zapis pozycyjny ze standaryzacją parametrów z wyd. 3.16 dla następującego zadania z trzema parametrami (podane maksymalne i minimalne wartości parametrów oraz długości podciągów kodowych):
max,: 20, 100, 300
min,.: -10, -5, 0
/,: 5, 10, 15
Wykonaj testy dla ciągu złożonego z samych zer, ciągu złożonego z samych jedynek oraz ciągu 010101 ... 0101. Porównaj wyniki komputerowe z obliczeniami odręcznymi. Określ także dokładność każdego z bloków składowych. Po przetestowaniu włącz rozważane podprogramy do programu SGA, pamiętając o procedurach ini-cjalizacji.
C. Zminimalizuj funkcję/(*,y,z)=*2+y2+z2, gdzie x, y, z mogą zmieniać się w zakresie od -512 do 512. Zastosuj kod lO-bitowy dla każdego parametru ^est to funkcja F\ De Jonga).
D. Zaimplementuj algorytm poszukiwania binarnego w oparciu o skumulowane prawdopodobieństwa selekcji w celu zwiększenia efektywności procedury reprodukcji.
E. Zaimplementuj i przetestuj podprogram zegara mutacyjnego do przeprowadzania mutacji. (Wskazówka: Posłuż się rozkładem geometrycznym; generuj "czas" do następnej mutacji).
F. Napisz podprogram realizujący kodowanie zmiennopozycyjne o zadanej wielkości mantysy i cechy.
G-. Porównaj efektywność algorytmu genetycznego z kodem dwójkowym i niedwójko-wym. W szczególności, porównaj efektywność AG z kodem dwójkowym i z kodem ósemkowym dla funkcji przystosowania f(x)-x10. Zastosuj 30-pozycyjny zapis dwójkowy i lO-pozycyjny zapis ósemkowy. Porównaj i omów tempo zbieżności i końcowy poziom zbieżności dla obu metod kodowania.