S U]\ ]D E X U]R Q \ FK Z DUWRĞFLDFK IXQNFML

52&=1,., 32/6.,(*2 72:$5=<67:$ 0$7(0$7<&=1(*2

6HULD ,,, 0$7(0$7<.$ 67262:$1$ ;;;9,,, 

5R E H U W :L H F ] R U N R Z V N L

:DUV]DZD

$ OJ R U \ WP \  VWR F K D VW\ F ] Q H  Z  R S W\P DOL]DFML G \VN UHWQ HM

S U]\ ]D E X U]R Q \ FK  Z DUWRĞFLDFK  IXQNFML

3UDFD ZSá\QĊáD GR 5HGDNFML 

$ E VWU D F W 6 WR F K D VWLF  D OJ R U LWK P V LQ  G LVFU HWH R S WLP L] D WLR Q  Z LWK  Q R LV\  YDOX HV

R I IX Q F WLR Q 7KH SDSHU GHDOV ZLWK VWRFKDVWLF PHWKRGV IRUVHDUFKLQJ DSSUR[LPDWHO\ JOREDO

PLQLPXP RIIXQFWLRQ GHILQHG RQ GLVFUHWH VHW $PHDVXUH RITXDOLW\ RIVROXWLRQ LV GHILQHG WR

FRPSDUH GLIIHUHQW DOJRULWKPV 6LPSOH 0RQWH &DUOR PHWKRG LV DQDO\VHG DV PDLQ DOJRULWKP

IRUZKLFKIRUPXODV GHDOLQJ ZLWK WKHPHDVXUHRITXDOLW\ DUHGHULYHG WZRFDVHV H[DFW YDOXHV

DQG QRLV\ YDOXHV RIIXQFWLRQ 7KLV 0RQWH &DUOR PHWKRG LV XVHG DV D EDVH LQ VLPXODWLRQ

H[SHULPHQWV IRU FRPSDULQJ RWKHU VWRFKDVWLF DOJRULWKPV 7KH VHFRQG SDUW RI WKH SDSHU

DQDO\VHV DV\PSWRWLF SURSHUWLHV RIWKH JHQHUDOL]HG VLPXODWHG DQQHDOLQJ DOJRULWKPV 7KHRU\

RI0DUNRYFKDLQVLVXVHGLQPRGHOOLQJ WKLVFODVVRIDOJRULWKPV 7KHRUHPVDERXW FRQYHUJHQFH

RIWKH UHFRUGV RIDOJRULWKPV WRVHW RIRSWLPD ZLWK SUREDELOLW\ RQH DUH SUHVHQWHG LQ WKHFDVH

RIIXQFWLRQ KDYLQJ UDQGRP QRLV\ YDOXHV 7KH SDSHU DOVR UHYLHZV NQRZQ UHVXOWV LQ WKHILHOG

RIVLPXODWHG DQQHDOLQJ W\SH DOJRULWKPV IRU IXQFWLRQ ZLWK UDQGRPO\ SHUWXUEHG YDOXHV

 6IRUPXáRZDQLH SUREOHPX

 : SURZDG]HQLH

'DQD MHVW IXQNFMD    ;  ² 5  JG]LH ;  MHVW ]ELRUHP G\VNUHWQ\P 3U]\M

PLHP\ ĪH ;   ^  1 `

1LHFOL IL  I L  GOD L … ;  RUD] QLHFK

I P L Q   P L Q I L  

L…;

1LHFK ;  R]QDF]D ]ELyU PLQLPyZ JOREDOQ\FK IXQNFML   F]\OL ;   ^L 

IL  IPLQ` =DGDQLH SROHJD QD Z\]QDF]HQLX  P Q RUD] SU]\QDMPQLHM MHG

QHJR] SXQNWyZLPLQ … $  3UREOHP\ PLQLPDOL]DFML QD]ELRU]H G\VNUHWQ\P

]NWyU\PLVSRW\NDP\ VLĊZSUDNW\FH GRGDMąGRWDNRJyOQLH VIRUPXáRZDQHJR

]DGDQLDSHZQHRJUDQLF]HQLH DPLDQRZLFLHOLF]QRĞü 1 ]ELRUX $MHVW ] UHJXá\

EDUG]R GXĪą OLF]Eą L Z SUDNW\FH QLH MHVW PRĪOLZH ]UHDOL]RZDQLH QDU]XFDMą

cego się rozwiązania., polegającego na przeglądaniu wszystkich N wartości funkcji. Jako modelowy przykład rozważymy problem komiwojażera.

P r z y k ła d 1. Danych jest m miast i symetryczna macierz odległości między nimi D = (d.^). Trasę komiwojażera definiujemy jako zamkniętą drogę w której każde miasto jest odwiedzone dokładnie raz. Zadanie polega na znalezieniu trasy o najmniejszej długości. Każda trasa może hyc rozpa- trywana jako element zbioru wszystkich cyklicznych permutacji m elemen- tów. Wynika stąd, że zbiór X w tym problemie ma liczność N = -- - W praktycznych zadaniach interesują nas wartości m > 100 , a to oznacza, że N > ^4, zatem N jest liczbą co najmniej rzędu 4.66631 • 10155.

Przyjmujemy zatem dalej, że ustalona jest liczba n, dużo mniejsza od N , która określa górne ograniczenie na liczbę punktów przeglądanych przez konkretny algorytm minimalizacji. Oczywiste jest, że przy takim ogranicze- niu problemu wyjściowego uzyskać można tylko rozwiązania przybliżone.

Rozważmy pewien podzbiór {żj,..., in} C X . Niech

= . min. fi-

Przez i* będziemy oznaczali dowolny spośród punktów {*i,..., in] dla któ- rego zachodzi /* (« i,..., i„) = fi*n. Wartość przyjmujemy za przybliżenie nieznanej wartości fmin, natomiast ż* za przybliżenie imin.

Zadanie wyjściowe sprowadza się do zadania takiego wyboru podzbioru { ?!,..., in} przez konkretny algorytm, aby wartość /* (ii,..., in) była, od- powiednio „bliska” wartości optymalnej f min-

Jako miarę jakości przybliżenia przyjmujemy następującą wielkość An — An{%\, • ■ •, in) — # {i ę X : fj < f* ( i i,...,in)} N

Liczba An określa względną wielkość zbioru punktów, dla których funkcja / przyjmuje wartości mniejsze od przybliżenia /*(«i,... ,in) wartości mini- malnej f min. Powyższa definicja nawiązuje do definicji rozważanej w książce [ZIN 86 ] (por. także [ZIE65]) dla przypadku funkcji określonej na ograniczo- nym podzbiorze R k‘, miała, ona ogólniejszą postać

_ X{x : f{x) < / « ) }

gdzie A oznacza miarę Lebesgue’a w R k, a z* oszacowanie minimum uzy- skane po n krokach algorytmu.

Wielkość An ma następującą własność: jeśli f* = f min, to An — 0;

w przeciwnym razie An > 0.

W pracy ograniczymy się do algorytmów stochastycznych, które polegają

na tym, że zbiór {zj,.. .,i n} wybiera się w pewien losowy sposób. Dokład-

nie mówiąc oznacza to, że algorytm stochastyczny utożsamiamy z miarą

Algorytmy stochastyczne iv optymalizacji dyskretnej

probabilistyczną P określoną na przestrzeni zdarzeń elementarnych A'n = I

I x ... x X .

'---v---' n razy

Podejście takie uzasadniamy następującymi argumentami:

• nie marny żadnych dodatkowych informacji dotyczących funkcji, czego wymagają metody deterministyczne;

• liczne algorytmy stochastyczne, stosowane do zagadnień praktycznych (np. zadania inżynierskie), dawały zadowalające rezultaty i okazywały się konkurencyjne w stosunku do algorytmów deterministycznych;

• algorytmy stochastyczne są łatwe w implementacji i pozwalają na. sze- rokie użycie komputerów.

Dla algorytmów stochastycznych szukania, minimum globalnego funkcji podzbiory {z i ,..., ?'n} są. losowe, zatem miara precyzji An jest zmienną lo- sową, której rozkład prawdopodobieństwa zależy od metody losowania punk- tów i i ,..., in, tzn. od P.

Ustalmy dwie liczby 7 /, 7 £ (0,1) oraz niech, dla ustalonej funkcji g : X n - P \ '

1 (U ł • • • 5 *'n ) • ^{'7 £ A . f i < f/(Z 1 , . . . , 7*n)}

< n Zdefiniujmy następującą wielkość

P(An < 11) = P(/l,?(/,:)).

Dla ustalonego algorytmu sensowne jest postawienie problemu: wyznaczyć najmniejszą liczbę n = n(N, //, 7 ) taką, że zachodzi

P{ kn < i]) > 7 .

Powyższą nierówność interpretujemy następująco: po wykonaniu n kroków algorytmu minimum globalne funkcji / jest z prawdopodobieństwem co naj- mniej 7 zlokalizowane z dokładnością do zbioru o względnej mierze nie więk- szej niż ?/. Rozwiązanie tak postawionego zadania w ogólności nie jest moż- liwe. Postąpimy zatem w sposób następujący: znajdziemy rozwiązanie po- wyższego problemu dla pewnego szczególnego algorytmu (prostej metody Monte Carlo) i posłuży nam to do dalszych porównań (symulacyjnych) za- chowania. się innych algorytmów optymalizacji.

1.2. Funkcje różnowartościowa i nieróżnowartościowa

Załóżmy, że funkcja. / nie jest różnowartościowa i liiecłi / będzie wielo- krotną. wartością tej funkcji. Niech A C A" będzie zbiorem, na. którym funkcja / przyjmuje wartość /. Niech rn = m ax{/(r( ) : f ( x t) < /} ora.z M = min{/(:/;,;) : f { x t ) > /}. Zmodyfikujmy dowolnie funkcję / do funkcji różnowartościowej przyjmującej na zbiorze X_ wartości w przedziale (

M).

Powtórzmy tę operację dla wszystkich wielokrotnych wartości / funkcji /. Z

definicji wielkości A„ wynika, że dla nowej funkcji przyjmuje ona wartości nie mniejsze, niż dla wyjściowej funkcji nieróżnowartościowej. W konsekwencji jeśli wyznaczymy n spełniające warunek

F{ A« < n) > 7

dla funkcji różnowart.ościowej, to wartość ta będzie spełniała, powyższy wa- runek również dla wyjściowej funkcji nieróżnowartościowej.

Wobec tego wszędzie dalej w pracy założymy o funkcji /, że jest różno- wartościowa, co w szczególności oznacza, że — 1. Odnotujmy jednak, że niektóre twierdzenia dotyczące algorytmów' sekwencyjnych będą obowią- zywały bez tego ograniczenia.

1.3. Funkcja z losowymi zaburzeniami

Rozważać będziemy też przypadek, gdy wartości funkcji są obliczane z pew- nym błędem losowym. Oznacza to, że po wylosowaniu podzbioru {żj,. . ., /n}

jako oszacowanie minimum przyjmujemy fn = m in j/i,.

gdzie fj = fj + €j oraz ej oznacza losowe zaburzenie. Zakładamy ponadto, że zaburzenia c?- są niezależnymi zmiennymi losowymi o jednakowym roz- kładzie, skończonym drugim momencie, własności E(ej) = 0, z pewną dys- trybuantą Fa £ X,

T — | Fa : Fa(x) = G , er > 0 parametr skali | ,

gdzie G jest dystrybuantą ustalonego rozkładu np. normalnego lub jedno- stajnego.

1.4. Organizacja pracy

Dalsza część pracy podzielona jest na dwa główne rozdziały: rozdział 2 doty- czący przypadku bez zaburzeń oraz rozdział 3 rozważający przypadek funk- cji z losowymi zaburzeniami. W podrozdziałach 2.1 i 3.1 wyprowadzamy wzory dotyczące zdefiniowanej miary jakości algorytmu dla prostej metody Monte Carlo. Stanowi to podstawę do porównań z innymi algorytmami lo- sowego szukania minimum funkcji, opisanymi w podrozdziale 2.4.

W podrozdziałach 2.2 i 2.3 analizujemy klasę sekwencyjnych algoryt-

mów optymalizacji, podając twierdzenia dotyczące zbieżności algorytmów

oraz odpowiednio ciągów rekordów. Podrozdziały 3.2, 3.3, 3.4 stanowią roz-

winięcie tych zagadnień dla przypadku z zaburzeniami losowymi. Na końcu

rozdziałów 2 i 3 (podrozdziały 2.5 i 3.6) przytaczamy wyniki symulacji dla

wybranych funkcji testowych, porównując je z uzyskanymi wcześniej wy-

nikami teoretycznymi. Analiza prostej metody Monte Carlo jest punktem

odniesienia do porównań (za pomocą symulacji) efektywności innych algo-

rytmów stochastycznych.

Algorytmy stochastyczne w optymalizacji dyskretnej

2. Przypadek bez zaburzeń 2.1. Prosta metoda Monte Carlo

Przejdziemy teraz do szczegółowej analizy konkretnego algorytmu stocha- stycznego tzw. prostej metody Monte ( 'ario. Metoda ta polega na tym, że losujemy kolejno punkty t j , ..., in ze zbioru X stosując schemat bez zwra- cania, przy czym zakładamy, że wylosowanie każdego podzbioru { h ,. .. ,?'n}

jest jednakowo prawdopodobne i równe (^) . Tak prosty schemat loso- wania daje w tym przypadku możliwość wyprowadzenia wzorów na miarę dokładności przybliżenia oraz analizę zagadnienia wyznaczania ?7(

V, 77 , 7 ).

Stanowić to będzie podstawę do porównań z innymi, bardziej złożonymi algorytmami losowego szukania minimum funkcji.

(₁)

L

1. Dla prostej metody Monte i -

i

N - d - l

71

Carlo zachodzi następujący wzór - , dla N — d — 1 > n

dla N — d — 1 < n

gdzie d = [Nrj\ oraz [

J oznacza największą liczbę całkowitą mniejszą lub równą x.

D ow ód. W świetle założenia o różnowartościowości funkcji /, bez zmniejszania ogólności można założyć, że /j < f 2 <■■• < I

- Niech By = {1,.. ., d+ 1}. Zauważmy, że interesujące nas zdarzenie losowe Ar]( f *) ma postać

{{hi • • • 1 Li} : 3(fc :! < / .’ < n) i^ G Bd}.

Dla zdarzenia przeciwnego łatwo otrzymujemy

, dla N — d — 1 > ri dla N — d — 1 < ??,,

skąd uzyskujemy wzór ( 1 ). ■

P(Ay(fn)1 = < (£) { 0

Dla danych ^ , 77,7 potrafimy zatem dobrać n = n{ N, //, y ) tak, żeby lJ( A,i < //) > 7 .

Okazuje się, że przy ustalonych 7 i //, n rośnie bardzo wolno gdy N — oc, i można wskazać takie

Y 0 ( 7 , 7/), że dla wszystkich N > A" 0 ( 7 , //)

n

N

7 7 , 7 ) =

gdzie

7 , 77 ) = liniA /-_+00

(i\r

7 ) (istnienie granicy wynika z dalszych

rozważań).

Ilustruje to następująca tabelka

Wartości n( N, ?/, 7 )

7 ⁷⁷

103 104 105 106 107 108 N 00 0.90 0.010 188 225 229 230 230 230 230 0.001 684 1888 2253 2297 2301 2302 2302 0.95 0.010 238 291 298 298 299 299 299 0.001 777 2383 2922 2987 2994 2995 2995

(

2

Dla wartości n( 00 ,?/, 7 ) mamy następujący wzór ln(l - ą r ( 0 0 , 77, 7 ) =

ln( 1 - 77 ) Wyprowadzenie tego wzoru jest następujące:

Wyprowadzimy najpierw wzór przybliżony dla wielkości P(Ar?(/*)), opierając się na wzorze Stirlinga

kl = V2irk kk exp(—k + d), 0 < 0 < — J_.ZrC 1 Wprowadźmy jeszcze oznaczenia:

Pd = c - t ' ) a ’ L(k) = klnk — k + -ln(27T k).

Wtedy

P(Avu: ) ) = i - n , L(k) < h\(k\) < L(k) + ₁₂ A:

Rozważane prawdopodobieństwo P(Av(f*)) będziemy przybliżali przez wielkość:

o ) p a?r(Ari(f:)) = i - p d,

gdzie pd jest przybliżeniem prawdopodobieństwa pd, określonym wzorem pd = exp(ln(pd)j,

przy czym ln(pd) jest następującym przybliżeniem dla ln(prf):

ln(pd) - L(N - c l - 1 ) + L(N - n) + - ¹

1

12 (iV — cl — 1 ) 12 (iV — n

- L(N) - L(N - cl - 1 - n).

Algorytmy stochastyczne w optymalizacji dyskretnej

Zauważmy, że błąd bezwzględny przybliżenia ln(pd) przez ln(pd) jest rzędu -Ąj. Stąd mamy

O < Pd ~ Pd

Pd < exp ₂

N ^-1

⁰

Zatem błąd bezwzględny przybliżenia prawdopodobieństwa P{Av(f*)) przez (3) jest rzędu -Ąj. Jeśli pd < | to również błąd względny tego przy- bliżenia jest rzędu gdyż wtedy

p[ An( p ) ) - Pd ~ Pd

P ( M f n ) )

Pd

Zauważmy, że

czyli także M Pd) > lii(Pd) Pd > Pd oraz

P(Av(f:)) > P ^ A ^ f * ) ) . Jeżeli zatem wyznaczymy n takie, że

p apr(A„(f:)) > 7 to również

PiAyif,:.)) > 7-

Fakt ten możemy wykorzystać do obliczeń n(N, ?/, 7 ) w przypadku umiar- kowanie dużych zadań tzn. dla N < N

(^, ?/), gdy obliczenia, według wzorów dokładnych byłyby zbyt kłopotliwe.

W celu wyznaczenia wartości ? 7 ( 7 V, 77 , 7 ) dla których, przy dużych N,

(4) P apr(Av(f *)) > 7 ,

postąpimy w następujący sposób.

Obliczymy granicę

lim P apr(Av(f*)) = lim (1 - pd) = lim (1 - exp(ln(pd))).

A — ►CC* Ar — ► CO A T—* 0 0

Mamy

lim ln(pd)) = lim -{ ( 2 iY (77 — 1 ) -f 2n + 1 )ln( N{ 1 — 77 ) — n — 1 ) N—

) N—t’-sc 2

+ (2.V - 2/7, + 1) 1

(,Y -

- (2N + 1 )ln(iV) - (2N(i] — 1) + l)ln (# (l - 77) - 1)}

= 11 ln(l — 77 ).

lim Pa^r(Arj(f*)') — 1 (1

Y

Zatem

Podstawiając obliczoną granicę do nierówności (4) otrzymujemy dowo- dzony wzór ( 2 ).

2.2. Algorytmy sekwencyjne

Wprowadzimy teraz ogólny schemat dla szerokiej klasy sekwencyjnych al- gorytmów optymalizacji. Ciąg punktów generowanych przez dany algorytm będziemy traktować jako pewien ciąg (W, k > 1 ) realizacji zmiennych lo- sowych o wartościach w A', który zostanie dalej zdefiniowany. Rozważamy następujący schemat:

Na wstępie przyjmujemy jako dane

• ciąg (R k, k > 1 ) niezależnych zmiennych losowych o jednakowym roz- kładzie jednostajnym na odcinku [ 0 , 1 ]

• ciąg liczb nieujemnych (Tkl k > 1)

• funkcje a,/3 : R X R+ — [0,1]

• macierz G = (gtj), i,j G A', gdzie gtJ > 0 oraz V?' Y,j=\ 9ij = J- Za pomocą elementów macierzy G definiujemy „otoczenia” X , C A”

punktów i G X \

Viex A i = {j G A : gij > 0 ).

Niecli mi = tt'i'g min{/(.T;) : x G Afi}. Punkt nazwiemy minimum lokal- nym funkcji /. Wielkość g%1 jest prawdopodobieństwem wyboru jako kandy- data do następnego kroku punktu j , jeśli algorytm jest w punkcie i.

Podamy teraz konstrukcję ciągu zmiennych losowych (Yk,k > 1). Niech Ij będzie wybrane w sposób dowolny. Jeśli już mamy elementy I j , ..., Ifi., to następująco definiujemy zmienną losową Yk+i: według ustalonego rozkładu prawdopodobieństwa {gyk j , j G A'} losujemy j G A' i definiujemy

{ 3 gdy fj - fvk > o i R k < a(fj - f n . Tk ) (A) Yk+1 = l j gdy fj - f Yk < 0 i R k < (3(f) - fyk,Tk)

[ Yk w przeciwnym przypadku

O zmiennej losowej R k zakładamy dodatkowo, że jest niezależna od zmiennych 1 j , ,. • •, Yk.

Interpretacja tej konstrukcji jest następująca. Algorytm startuje z lo-

sowo wybranego punktu łfi zbioru A'. Przypuśćmy, że w k-tym kroku algo-

rytm jest w punkcie Yk = i. Wtedy z prawdopodobieństwem gtJ wybieramy

punkt j jako „kandydata” na kolejne przybliżenie Yk+\. Nie zawsze jednak

tego kandydata akceptujemy (jeżeli nie ma akceptacji, to Yk+ \ = Yk ). Ak-

ceptacja ma charakter losowy. .Jeżeli kandydat kandydat reprezentuje lepsze

rozwiązanie, tzn jeżeli Yt < f Yk, akceptujemy go z prawdopodobieństwem

(3 = (3(fj — fyk, Tk); jeżeli reprezentuje gorsze rozwiązanie, akceptujemy go

z prawdopodobieństwem a = cĄj) — fyk, Tk). Ta wstrzemięźliwość przy ak-

ceptacji lepszego rozwiązania lub cofanie się do rozwiązania gorszego ma,

uchronić algorytm przed „ugrzęźnięciem” w punkcie minimum lokalnego.

Algorytmy stochastyczne w optymalizacji dyskretnej

Możemy teraz policzyć następujące prawdopodobieństwo (5) Pr{Yfc+i = j\Y[ = 2 / 1 , , Yk-i = Uk-i, Yk = i}

__ f gijFu{a(fj - fi,Tk)}, jeśli f, - f t > 0 1 gijFu{(3{ fj - fi, Tk)}, jeśli fj - f t < 0 gdzie Fu jest dystrybua.ntą rozkładu jednostajnego na [ 0 , 1 ].

Z konstrukcji zmiennych losowych Yk i wzoru o wynika, że ciąg {Yk,k >

1 ) jest łańcuchem Markowa z przestrzenią stanów X i macierzą przejść w jed- nym kroku P(k) = [Pij{Tk)]iJeX, gdzie

P.l:j{T) = PrjYfc+i = j\Yk = ?.} = gij(itj{T), śli fj - fi > 0

fi,T)}, jeśli fj - fi < 0

„ i T, \ F u { a U j - f i , T ) } ,

ij[ } ~ i F u m / j Yi.j e x , r > o .

Wielkość (iij{Tk) określa się jako prawdopodobieństwo akceptacji punktu j, jeśli w k-tyin kroku znajdowaliśmy się w punkcie i. Funkcję «^-(T) bę- dziemy dalej określali terminem funkcja akceptacji. Odnotujmy, że zależy ona od pewnego parametru T, którego wartości w kolejnych krokach algo- rytmu są wyznaczone przez ciąg (Tk, k > 1). Parametr T (a także ciąg Tk) jest więc parametrem „sterującym”. W algorytmach typu „simulated anne-

aling” interpretowany jest jako „temperatura” pewnego procesu fizycznego;

mówimy o tym dokładniej w rozdziale .

Potrzebne nam będą jeszcze następujące oznaczenia:

a(P) = min i, l L ' j mh\(Pi j, Pij)

(jest to tzw. współczynnik ergodyczności dla macierzy stochastycznej P, patrz np. Iosifescu 1988, rozdział 1 . 11 , str. 54).

Niech M będzie taką najmniejszą liczbą naturalną, że dowolny punkt należący do zbioru punktów optymalnych A'* może być osiągnięty z każdego innego punktu po co najwyżej M krokach, czyli

M = min {Z : Vi€A-,je _

* fJ^j > °K gdzie gf - jest elementem macierzy

GL = G x G x .. .x G .

L razy

Wprowadzone wyżej elementy (Rk,k > 1 ), (Tk,k > 1 ).

, fj oraz G

w pełni wyznaczają łańcuch Markowa {Yk,k > 1 ). Zamiarem konstrukcji

jest, aby łańcuch był zbieżny ( w sensie, o którym dalej powiemy) do zbioru

minimów globalnych A'*. Może sie jednak zdarzyć, że ten łańcuch ma rów-

nież inne stany powracające (w sensie np. definicji Iosifescu 1988, rozdział

2.5). Stany powracające, które są jednocześnie minimami lokalnymi, a nie

są minimami globalnymi, będziemy nazywali punktami suboptymalnymi, a.

zbiór tych punktów oznaczymy przez A'*.

Niech

a(k) = min( min ciij(Tk)) a(k) = max( max fi i, {'/),)).

~ iex je-

iex» jeXi\x*

Niech Tl{T) oznacza rozkład stacjonarny związany z macierzą Wprowadźmy ponadto oznaczenie

p ( m , k )

P ( m ) .. .P(k) jeśli m < k I jeśli m > k

W pracy [ANLSTaJ udowodniono twierdzenie dotyczące warunków ko- niecznych i dostatecznych na to aby algorytm (A) miał następujące włas- ności

(a,) Osiągalność zbioru optimów globalnych A* (mówimy, że zbiór optimów globalnych A'* jest osiągalny, jeżeli dla każdego Ij ciąg l j , y 2, . .. zbiega z prawdopodobieństwem 1 do A'*)

(b) Asymptotyczna niezależność od rozwiązania początkowego, czyli tzw. słaba ergodyczność łańcucha Markowa oznaczająca,że

Vi,j,leX, m> 1 lim

k —+oo

p{ m,k)

1 ij

T)(m,k)

1 lj

=

(c) Zbieżność Yk według rozkładu, równoważna tzw. mocnej ergodyczno- ści łańcucha: istnieje graniczny rozkład II = (III, IN ,. .., II/v) taki, że

V , m > ! lim

—

o P p k) = n, (d) Zbieżność z prawdopodobieństwem 1 do zbioru A'*.

Twierdzenie to zakłada spełnianie przez funkcje akceptacji (iij{T) nastę- pujących warunków regularności:

(iij{T) \ Oclla T dostatecznie małego, jeśli fj > fi ]__ ciij[T) Idla T dostatecznie małego, jeśli fj < fi hm ciij (T )istnieje, jeśli fi = f r

2. Wektor prawdopodobieństw rozkładu stacjonarnego 11(7’) jest •funkcją (wektorową) R + 3 T U(T) 6 RN o wahaniu ograniczonym.

Funkcjami spełniającymi powyższe warunki są np. (patrz tw. 2 ,

[ANI87b]) funkcje wymierne wielomianów lub funkcje wykładnicze zmien-

nych T lub T -1 , lub funkcje będące sklejeniem funkcji należących na odcin-

kach do powyższych klas.

Algorytmy stochastyczne w optymalizacji dyskretnej

T

1 . i) Jeśli CLM{kM) = oo, to zachodzą własności (a),

(b) i (c) oraz istnieje lirii/i:_,Xl FI (k). Jeśli ponadto lim n(^)j = 0 Vj;/j>/l

k ---T OO

to zachodzi również warunek (cl).

ii) Jeśli istnieją punkty suboptymalne i '52‘kLi dM{k) < oo to nie zachodzi żaden z warunków (a)-(d).

Dowód lego twierdzenia znajduje się w pracy [ANI87a]. ■ Przykłady funkcji akceptacji spełniających warunki regularności:

1) funkcja wykładnicza, (Metropolisa): af'jet(T) = min j 1, exp( ) j . 2 ) funkcja Hastings a:

1 + 210.5 min{ ~ exp( L ^ 9j i fi -Ii

T

9 ii

9ij ^exp( fi-fi

)}]'

1 + ^ exp( 9j gdzie 7 > 1 jest ustalonym parametrem.

Dla obu rodzajów funkcji znane są wyrażenia dla IF(Tfc) i wiadomo, że lim ^oo 11(/r), = 0 , gdy j nie jest globalnym minimum ([HAS.70] i [LT.JN 86 ]).

Niech A+ = max{ fj — ft : j G Xj,i 6 A'} oraz A = min{ /j — ft : i ę x \ j e X i \ x * } .

Twierdzenie 1 zastosowane do powyższych funkcji implikuje, że warunki (a)-(cl) zachodzą, jeśli J\ > dla k — oo, podczas gdy żadna z roz- ważanych własności nie zachodzi dla żadnego problemu z punktami subop- tymalnymi, jeśli tylko Tk < , A: -* oo. Oznacza to, że dla obu typów funkcji akceptacji zachodzenie warunków (a)-(d) ma miejsce, jeśli zmiana parametru Tk następuje dla dostatecznie dużych k nie szybciej niż

Dodatkowo z twierdzenia 1 wynika cieką,wy wniosek dotyczący tego, że dla wielu popularnych schematów zmiany parametru Tk, np.

Tk = /)k, ji < 1,

nie zachodzi żadna z własności (a)-(d), gdy stosujemy ten schemat do pro- blemu z punktami suboptymalnymi. Odmiana algorytmu w której Tk jest zmieniane po określonej liczbie iteracji {Aą}^ (tzn. Tkl — ckl+ 1 = ... = Tkl+1-i) nie spełnia (a)-(d), jeśli

Tk i = I~e oraz Aą = V/ 3 >0,e> o gdyż wtedy

CO OO

^ d ( T r ) = ^ / ^ e x p ( - A 'r ) <

v=l i= i oo.

2.3. Ciągi rekordów

Oprócz ciągu (Yk, k > 1) interesować nas również beclzie ciąg tzw. rekordów algorytmu (A), określony w sposób następujący:

rk - arg min f(Ym), k = 1 , 2 ,...

1 < m < k

Rozważanie ciągu {rk,k > 1) można uzasadnić tym, że w praktycznych realizacjach algorytmu musimy zawsze kiedyś zakończyć obliczenia i wtedy sensownym przybliżeniem szukanego minimum może być właśnie wartość rekordowa wśród dotychczasowych wartości ciągu Yk. Naszym celem będzie pokazanie przy spełnieniu jakich warunków ciąg rekordów zbiega z prawdo- podobieństwem 1 do zbioru A’*.

Ponieważ rozważany zbiór A” jest zbiorem dyskretnym, zachodzi nastę- pujący oczywisty fakt:

L^{e m a t} 2.

Przy spełnieniu warunków twierdzenia 1 zachodzi zbieżność

ciągu rekordów (rk,k > 1 ) z prawdopodobieństwem 1 do zbioru X *.

2.4. Algorytmy optymalizacji globalnej

Będziemy rozważali następujące algorytmy optymalizacji globalnej:

2.4.1. Algorytm MC. Przez algorytm MC rozumiemy algorytm prostej metody Monte Carlo opisany w rozdziale 2.1.

2.4.2. Algorytmy SA. Algorytmy SA oznaczały będą klasę algoryt- mów typu symulowanego wygrzewania, (simulated annealing) (patrz m.in.

[ANI87a, BÓN84, CER85, JMM 88 , KIR.82, KIR83, LAA87]). Zostały one skonstruowane w oparciu o idee zaczerpnięte z fizyki statystycznej. Kirkpa- trick [KIR.83] i Cerny [CER85] zauważyli pewną analogię między minimaliza- cją funkcji (zadaną bez błędów losowych) w problemach kombinatorycznych a powolnym ochładzaniem i utwardzaniem (wygrzewaniem) ciała stałego, aż osiągnie ono stan o minimalnej energii. Ponadto stwierdzono, że proces optymalizacji może być zrealizowany z wykorzystaniem metody Metropolisa.

[MET53] do symulacji ewolucji ciała, stałego w stanie równowagi termody- namicznej. Dzięki rozpowszechnieniu nowego pomysłu nastąpił istotny roz- wój badań teoretycznych nad nową klasą algorytmów oraz w dziedzinie jej zastosowań. W interpretacji algorytmu typu SA funkcja, którą chcemy zmi- nimalizować odpowiada, energii ciała stałego. W każdym kroku algorytmu system jest losowo przekształcany z aktualnego stanu do nowego, znajdu- jącego się w otoczeniu stanu aktualnego. Oblicza się zmianę w funkcji SE.

Jeśli następuje redukcja wartości funkcji, to przejście do nowego stanu jest

akceptowane, (we wzorze (A) /3 = 1), w przeciwnym przypadku akceptacja,

następuje tylko z pewnym prawdopodobieństwem, zwykle postaci exp( — ),

gdzie T jest parametrem kontrolnym, odpowiednikiem temperatury. Powyż-

szy schemat, w mechanice statystycznej nazywany schematem akceptacji

Algorytmy stochastyczne w optymalizacji, dyskretnej

Metropolisa, symuluje zachowanie się ciała stałego osiągającego równowagę termodynamiczną. Po odpowiednio dużej liczbie kroków rozważany układ w temperaturze T zbliża się do stanu równowagi tzn. rozkład stanów energii jest bliski rozkładowi Boltzmanna, czyli:

Pr {energia układu < E } = f

— CO

K f ) exp — x kBT ^{d x}

gdzie Z{T) jest czynnikiem normalizującym, a kp tzw. stalą Boltzmanna.

Z doświadczeń fizycznych wiadomo, że przy zbyt szybkim ochładzaniu (nie osiągana jest równowaga termodynamiczna w każdej temperaturze ) ciało nie uzyskuje struktury krystalicznej o minimalnej energii, natomiast utrwalają się w nim defekty odpowiadające lokalnym minimom energii. Problem szyb- kości zmniejszania parametru kontrolnego (temperatury) okazuje się rów- nież kluczowy, gdy przeniesiemy powyższy schemat do zadań optymalizacji kombinatorycznej: należy tak zmieniać parametr T, aby algorytm uniknął wpadnięcia w pułapki lokalnych minimów.

Poniżej przedstawiamy zapis algorytmu SA w pseudokodzie w postaci procedury anneal,

procedure anneal;

{używane funkcje i parametry:

accept(Ćif,T) - akceptacja lub odrzucenie nowego stanu T(t) - parametr temperatury z kolejnym numerem t N(t) - liczba powtórzeń przy ustalonej temperaturze T(t) } begin

inicjalizacja stanu i;

t:=0; { inicjalizacja licznika zmian temperatury } oblicz T(0); { inicjalizacja temperatury }

do while (kryterium stopu = false);

n: = l; { inicjalizacja licznika powtórzeń } repeat

generacja stanu j z sąsiedztwa i;

oblicz SE; {zmiana w funkcji celu } wywołaj funkcję accept( CE,T(t));

if (accept = true) then

i:=j; { j staje się stanem bieżącym } n: = n + l

until ( n = N(t));

oblicz T(t-fl);

t:=t + l;

endwhile;

end;

function accept(ÓjF ,T);

{używane funkcje:

random(0,l) - realizacja zmiennej losowej o rozkładzie jednostajnym z przedziału (0,1) }

begin

if (SE < 0) then accept:=true

else if (random(0,l) j exp( — 4^r) ) then accept:=true else accept:=false;

Konkretne implementacje powyższego schematu wymagają skonkretyzo- end;

wania tzw. planu ochładzania, na co składają się:

• wartość początkowa parametru T: duża. wartość na początku oznacza., że nowy kandydat, na kolejne przybliżenie jest przeważnie akcepto- wany, a. obniżanie T oznacza mniejsze prawdopodobieństwa, akceptacji gorszy cliś tanów

• funkcja zmiany temperatury T (t), t oznacza licznik zmian temperatury

• liczba powtórzeń w danej temperaturze N(t )

• kryterium stopu.

W literaturze można spotkać wiele wariantów planu ochładzania. Oto przykładowe możliwości:

Funkcja temperatury:

• stała: T(t) = C

• arytmetyczna: T(t) = T(t — 1 ) — C

• geometryczna: T(t) = a(t)T(t — 1), a(t) często jest stalą

• odwrotna: T(t) =

• logarytmiczna: T{t) = [n(f^T) Liczba, powtórzeń:

• pojedyńcza: N(t) = 1

• stała: N(t) = C

• arytmetyczna,: N(t) = N (t — 1 ) + C

• geometryczna: N(t) = —

• logarytmiczna: N{t) =

• potęgowa: N(t) = N(t — l)«óT

• energia: powtarzaj aż oszacowanie średniej energii układu •mało się zmieniło

• akceptacja: powtarzaj aż do określonej liczby akceptacji

• odrzucenie: powtarzaj aż do określonej liczby odrzuceń Kryterium stopu:

• iteracje: określona liczba iteracji

• temperatura: T(t) < ustalona warto kocowa T

Algorytmy stochastyczne w optymalizacji dyskretnej

• energia: oszacowanie średniej energii mało się zmienia w kolejnych tem- peraturach

• akceptacje: w kolejnej temperaturze występuje zbyt mala liczba ak- ceptacji

• poprawa: w kolejnych temperaturach nie było poprawy funkcji celu.

W niniejszej pracy w eksperymentach symulacyjnych zastosowano nastę- pujące parametry algorytmu SA:

• T{ 0) = 0.5,

• T(t) = 0.9 T[t — 1),

• N{t) = 10,

• kryterium stopu: ustalona liczba iteracji. Liczba ta jest postaci n =

??,(?/, 7 ) i wynika z założonej dokładności rozwiązania ( w sensie defi- nicji z rozdziału 1 . 2 ) dla. prostej metody Monte Carlo. Algorytm MC dzięki swej prostocie pozwala na łatwe obliczanie ??,(?/,. 7 ) i stanowi jednocześnie przykład modelowy, do którego odnosimy zachowanie się

innych algorytmów optymalizacji stochastycznej.

2.4.3. Algorytmy TA. W pracy [DUE90] zaproponowano pewną mo- dyfikację algorytmu SA, nazwaną „threshold accepting” (TA). Nowa klasa, algorytmów charakteryzuje się prostszą strukturą, niż SA; w rozważanych problemach komiwojażera i konstrukcji kodów binarnych uzyskiwano dzięki niej rozwiązania bliskie optymalnym w czasie o wiele mniejszym, niż dla algorytmów SA.

Jako parametr kontrolny T w algorytmie typu TA wykorzystywane są tzw. poziomy akceptacji (thresholds). Stanowią one po prostu ustalony nie- rosnący ciąg nieujemnych liczb rzeczywistych (Jj., k > 1 ). W realizacjach algorytmu poziomy wprowadza się jako skończoną tablicę liczbową. Sche- matycznie algorytm TA można zapisać w postaci następującej procedury TA:' procedure TA;

wybierz punkt początkowy;

nr = 1; { początkowy numer poziomów akceptacji } wybierz początkową wartość poziomu T := Tnr ; Opt: wybierz nowy punkt z otoczenia bieżącego;

oblicz ÓE = f(nowy punkt) - f(stary punkt) ; IF {ÓE < T) THEN stary punkt:= nowy punkt;

IF ( zbyt dużo iteracji dla danego poziomu T )

THEN ( nr:=nr+l; T := Tnr ) ; { obniżenie poziomu ) IF (spełnione kryterium stopu) THEN stop ELSE GOTO Opt;

Kryterium stopu można, definiować na wiele sposobów, podobnie jak dla

algorytmu SA. Najczęściej jako warunek zakończenia obliczeń przyjmuje się

brak poprawy dla ustalonej liczby iteracji. W niniejszej pracy w dalszych eks- perymentach symulacyjnych jako kryterium stopu przyjęto ustaloną liczbę iteracji (patrz wyżej).

Istotna różnica między schematami SA i TA polega na odmiennych re- gułach akceptacji nowego punktu o większej wartości funkcji celu. Algorytm TA akceptuje każdy nowy punkt, w którym wartość funkcji jest odpowiednio większa od starej, natomiast SA akceptuje punkt o większej wartości funk- cji tylko z pewnym prawdopodobieństwem. Zaletą TA jest jego prostota; nie wymaga się w nim obliczania prawdopodobieństw i podejmowania, losowych decyzji.

2.5. Porównania symulacyjne

Poniższy rozdział zawiera wyniki eksperymentów symulacyjnych przeprowa- dzonych dla konkretnych funkcji testowych i wybranych algorytmów opty- malizacji.

Ogólna idea eksperymentu była następująca. Dla ustalonych wartości i] oraz 7 , odpowiadających założonej dokładności rozwiązania, obliczamy liczbę n, = 77 ,( 77 , 7 ) według wzoru (2) str. 124. Jest to ilość punktów, jaką należy wylosować ze zbioru A" w schemacie prostej metody Monte Carlo, aby rozwiązanie miało zadaną dokładność. W przeprowadzonych ekspery- mentach symulacyjnych przyjęto // = 0.01 oraz 7 0.9, co oznacza, że minimum globalne funkcji ma być z prawdopodobieństwem co najmniej 0.9 zlokalizowane z dokładnością do zbioru o mierze względnej nie większej niż 0.01. Każdym z porównywanych algorytmów MC, SA oraz TA wykonujemy tą samą. ilość n „kroków iteracyjnych”, przy czym rozważamy po dwie wer- sje algorytmów SA i TA, różniące się między sobą wyborem macierzy G: we wszystkich wersjach jest to taka macierz, że otoczeniem Ab punktu i G A” jest

{(i — D ) mod N , ..., (-i + D ) mod A'},

gdzie D jest parametrem określającym wielkość otoczenia, a rozkład < 7 .;. na tym otoczeniu jest rozkładem jednostajnym. W wersjach SAa i TAa przyj- mujemy D = 102, a w wersjach SAb i TAI) D — 105. Oznacza to, że algorytm SAa (oraz TAa) ma, „bardziej lokalny” charakter niż algorytm SAb (odpo- wiednio TAb).

Po wykonaniu n kroków każdym z algorytmów rejestrowano fakt zajścia lub niezajścia zdarzenia Z\ „minimum globalne funkcji zastało zlokalizo- wane z dokładnością do zbioru o mierz: ..jględnej < //”. Każdy z takich eksperymentów został powtórzony M = 1U000 razy. Niech M(MG) będzie ilością tych eksperymentów z algorytmem MC, w których zaobserwowano zdarzenie Z. Analogicznie definiujemy M(SAa),M(SAb),M(TAa),M(TAb).

Wtedy A-(a-lgerytinJ jest empirycznym oszacowaniem prawdopodobieństwa

P(An < rj) (por. str. 121 ), i w przypadku algorytmu MC jest estymato-

rem znanej (założonej) liczby 7 = 0.9.

Algorytmy stochastyczne w optymalizacji dyskretnej

W podanych niżej tabelkach przedstawiono wartości 'V^ al^°/I>tm^ dla wszystkich pięciu rozważanych algorytmów, dla, dwóch funkcji testowych J

-, I

'■ di — Ii. Funkcja J

jest zdehniowana wzorem

J

{i) = i, i natomiast

f B{i) =

(/). ż = 1,.. .,iV,

gdzie 7r() jest pewną losową permutacją zbioru {1,2,..., N}.

N — 2 16, funkcja J

11 MC SAa. SAb TAa TAb

230 .90 .10 1.00 .10 1 . 00 - N = 218, funkcja J

n MC SAa SAb TAa TAb 230 .90 .03 1.00 .03 1.00

N = 216, funkcja }B n MC SAa SAb TAa TAb

230 .91 .77 .91 .78 .91

N = 218, funkcja f B n MC SAa SAb TAa TAb

230 .90 .77 .90 .78 .91

Z analizy eksperymentów symulacyjnych wynikają następujące wnioski:

Oszacowania uzyskane z symulacji dla prostej metody Monte Carlo potwier- dzają wyniki teoretyczne uzyskane dla tego algorytmu. Zbliżone rezultaty dla funkcji J

i f B wskazują też na małą wrażliwość metody Monte Carlo na zmiany postaci funkcji, co w naturalny sposób wynika z „globalnego"

mechanizmu działania tego algorytmu.

Nie zauważa się istotnych różnic w działaniu algorytmów typu SA i TA.

Dla efektywności obu algorytmów widoczny jest decydujący wpływ struk-

tury otoczeń (warianty a i b). Losowania typu globalnego (b) daje wyniki

zbliżone do uzyskanych przez metodę MC. Losowanie z otoczeń, które są

różne od całego zbioru A' powoduje bardziej „lokalne” działanie algoryt-

mów i stąd istotne zmniejszenie efektywności algorytmów SA i TA dla obu

funkcji testowych. W przypadku (a) przeprowadzony eksperyment potwier-

dza znaną z licznych doświadczeń obliczeniowych własność algorytmów typu

symulowanego wygrzewania: wymagają one bardzo długiego czasu obliczeń

(czyli dużego n), a jak się okazuje, w pewnych zadaniach dużo większego,

niż w algorytmie MC.

Podsumowując przeprowadzone eksperymenty obliczeniowe można stwierdzić, że mając, na uwadze efektywność algorytmu optymalizacji mie- rzoną zdefiniowaną miarą precyzji, nie uzyskano potwierdzenia dla rozpo- wszechnionego poglądu o zdecydowanej przewadze metod typu symulowa- nego wygrzewania.

Wydaje się również, że zdefiniowana miara, precyzji stanowić może dobre narzędzie porównawcze w dalszych badaniach różnych algorytmów optyma- lizacji dyskretnej. W szczególności istnieje możliwość jej użycia w badaniach symulacyjnych, w sytuacjach gdy trudno jest uzyskać wzory teoretyczne do- tyczące danego algorytmu.

WT powyższych badaniach użyliśmy funkcje testową fsi'i) = 7r(/), gdzie 7 r było „pewną losową permutacją” zbioru {1,2...., N }. Wyjaśnimy teraz szczegółowo konstrukcję tej permutacji.

Ponieważ interesujące są, zadania z dużymi wartościami N i przechowy- wanie w pamięci całej informacji o permutacji

) staje się niemożliwe, ko- nieczna, jest umiejętność algorytmicznego wyznaczania, wartości funkcji fs- W celu ułatwienia, obliczeń dla, funkcji fg przyjęto następującą konstrukcję:

wykorzystano generatory liczb pseudolosowych typu ( 6 )

— (a x-i -1 + 6 )mod 2 m, przyjmując

/

('0 = + 1 ,

gdzie Xi oznacza, ż-tą wartość z generatora,, który startuje z ustalonej warto- ści x0. Generator typu ( 6 ) posiada następujące własności (patrz up. [ZIE79], rozdz. 2, str. 47): maksymalny okres tego generatora, jest równy 2m, nie za- leży on od x'o oraz jest osiągalny dla a i 6 spełniających warunki

a — 1 mod 4, 6 = 1 mod 2.

Spełnienie powyższego warunku oznacza,, że ciąg (

l

2 m ) daje

róż- nych liczb ze zbioru {0,.. .,2 m — 1}. Zatem przyjmując w definicji funkcji J

odpowiednie parametry a, 6 ,ustalając x{) oraz N = 2 m uzyskujemy moż- liwość obliczania, wartości funkcji /s(i) „na bieżąco”, bez zapamiętywania, informacji o całej permutacji, czyli tablicy liczb ( /

(1), . .., /jg(lV)). Do- datkowe przyspieszenie obliczeń uzyskuje się tworząc na, początku tablicę przechowującą wybrane wartości Xi z generatora.

Do obliczeń użyto następujących układów współczynników:

m N = T n a 6 x0

16 65536 10001 77 100

18 262144 441 3333 10000

Algorytmy stochastyczne w optymalizacji dyskretnej

.Jako wartość startową dla badanych algorytmów przyjmowano punkt, wylosowany w sposób jednostajny ze zbioru X . Do generowania liczb pseu- dolosowych użyto generatora opisanego w pracy ([ZAM90]).

Dla przyjętej wielkości M = 10000 stwierdzamy, że z prawdopodobień- stwem co najmniej 0.9 błąd uzyskanych oszacowań nie przekracza 0.01.

2.6. Problem komiwojażera

Wśród wielu problemów optymalizacji dyskretnej ważnym zagadnieniem praktycznym jest wspomniany już na początku pracy problem komiwoja- żera. Duża złożoność obliczeniowa tego problemu (należy on do klasy tzw.

zadań NP - zupełnych) spowodowała, powstanie wielu szczególnych przybli- żonych metod jego rozwiązywania. Stanowi on ponadto przykład modelowy na którym weryfikować można przydatność nowych algorytmów optymali- zacji dyskretnej.

Przegląd klasycznych metod rozwiązywania problemu komiwojażera wraz z obszernymi omówieniami algorytmów i wyników doświadczeń ob- liczeniowych zawiera książka [SYS93].

Algorytm symulowanego wygrzewania stosowany do przykładowych pro- blemów komiwojażera dawał zadowala jące wyniki w porównaniu z metodami znanymi wcześniej. Szczegółowe omówienie doświadczeń i odniesienia do li- teratury znaleźć można w książce [LAA87].

Dla specyficznych problemów optymalizacji ważny jest mechanizm wy- boru elementów z sąsiedztwa danego punktu (czyli trasy). Trasę komiwoja- żera wygodnie jest traktować jako tablicę ROUTE[l..m+ 1 ], gdzie element ROUTE[i],i = 1 ,..., m j- 1 oznacza kolejny numer odwiedzanego miasta oraz ROUTE[rn + 1] := ROUT E[\], W algorytmach dla zadań komiwo- jażera stosuje się z reguły mechanizm przejścia wprowadzony przez Lina i Kernighana, są to tzw. przejścia k-opt. Najprostsza metoda 2 -opt polega na tym, że wybierane sa 2 miasta i oraz j z bieżącej trasy, a następnie zmie- nia się kolejność odwiedzania miast na odcinkach między miastami i oraz j. W ten sposób można, zdefiniować otoczenie jako zbiór takich tras, które mogą być osiągnięte z bieżącej poprzez przejścia typu k-opt.

Poniżej przedstawia,my schemat, przejścia, 2 — opt:

procedure 2-opt:

wybierz dwa miasta i,j

{1 ,.. . ,m }, i < j;

likwidujemy następujące połączenia:

między miastami ROUTE[i] i ROUTE[i + 1]

oraz między ROUTE[j] i ROUTE\j + 1]

łączymy miasta:

1 lOUTE[i\ z ROUTE[j]

oraz ROUTE[i + 1] z ROUTE[j + 1]

Standardowy algorytm TA dla zadania komiwojażera stosowany w pracy [DUE90] miał następującą postać (procedura, TAI):

procedure TAI;

wybierz losowo początkową trasę;

inicjalizacja tablicy poziomów

T[]:= (0.13,0.12,0.11,0.10,0.095,0.09,0.085,0.08,0.075,0.075, 0.075,0.07,0.07,0.07,0.065,0.065,0.065,0.065,0.06,0.06, 0.055,0.055,0.05,0.05,0.05,0.04,0.04,0.03,0.02,0) ; FOR każdego poziomu T=0.13,0.12,...,0

DO s razy { s - ustalona liczba powtórzeń } wybierz miasta i, j losowo;

2-opt; wybór nowej trasy

IF długość nowej trasy j długość starej + T THEN stara trasa:= nowa trasa ;

Zbadano również odmianę deterministyczną powyższego algorytmu (pro- cedura TA2):

procedure TA2;

ustal liczbę całkowitą b (w pracy b=10);

dla każdego miasta znajdź b najbliższych sąsiadów;

wybierz losowo trasę początkową ;

T[]:= (0.13,0.12,0.11,0.10,0.095,0.09,0.085,0.08,0.075,0.075, 0.075,0.07,0.07,0.07,0.065,0.065,0.065,0.065,0.06,0.06, 0.055,0.055,0.05,0.05,0.05,0.04,0.04,0.03,0.02,0) ; FOR każdego poziomu T=0.13,0.12,...,0

FOR każdego miasta i: = l,...,m

FOR każdego miasta będącego jednym z b najbliższych dla i

DO wykonaj przejście typu 2-opt ;

IF długość nowej trasy j długość starej + T THEN stara trasa:= nowa trasa ;

W deterministycznej wersji algorytmu TA w procedurze przejścia do no- wej konfiguracji metodą 2 -opt nie wybiera się 2 miast losowo, tylko w okre- ślonym z góry porządku.

Oba schematy zastosowano w pracy [DUE90] do znanego w literatu- rze problemu komiwojażera dla 442 miast z pracy [GR084], gdzie ana- lizowano rezultaty stosowania SA. Eksperymenty z algorytmem SA po- twierdzają, że kluczem do osiągnięcia dobrych rezultatów jest odpowiedni

„plan ochładzania”. Natomiast eksperymenty z TA wykazują, że jest on

raczej odporny na wybór poziomów akceptacji. Wniosek główny jest na-

stępujący: algorytm TA przewyższa jakością SA, gdyż średnio daje dużo

lepsze rozwiązania, niż otrzymane do tej pory przez algorytm SA. Dodat-

kowym zaskakującym wnioskiem wynikającym z eksperymentów jest fakt,

Algorytmy stochastyczne w optymalizacji dyskretnej

że deterministyczna wersja algorytmu TA daje równie dobre, bliskie opty- malnym rozwiązania w znacznie szybszym czasie. Czas wykonania standar- dowej wersji TA dla zadania z 442 miastami wynosi średnio kilka minut (około 4,000,000 iteracji), co jest porównywalne z wynikami dla SA, na- tomiast dla wersji deterministycznej czas ten skracał się do kilku sekund (442,000 iteracji); obliczeń dokonywano na komputerze IBM 3090 Mod 200 V F. Spośród konstruowanych w ostatnim czasie specyficznych algorytmów rozwiązania problemu komiwojażera z symetryczną macierzą odległości bar- dzo efektywnym okazał się algorytm T. Volgenanta i W.B. van den 11 out a z Institute of Actuarial Science and Econometrics, University of Amster- dam. Metoda zasosowana w nowym algorytmie opiera się na pojęciu 1-drzew i wykorzystuje schemat podziału i ograniczeń oraz relaksację lagrangowską ([VOL82]).

W niniejszej pracy przeprowadzono dodatkowo następujący ekspery- ment: dla losowo wybranego układu miast reprezentowanych przez m punk- tów na płaszczyźnie w kwadracie jednostkowym wykonano n iteracji al- gorytmami prostej metody Monte Carlo oraz algorytmem SA, przy czym uwzględniano dwa mechanizmy wyboru sąsiedztwa dla bieżącej trasy komi- wojażera w algorytmie SA

• RANPER. - polegający na losowym wyborze permutacji ciągu ( 1 ,..., ???.), przy czym prawdopodobieństwo uzyskania danej permu- tacji równe jest ^ ,

• 20PT - opisany wcześniej mechanizm typu 2-opt.

Trasa początkowa we wszystkich przypadkach była. ustalana w sposób losowy. Obliczenia powtarzano 10 razy. Poniższa tabela przedstawia wyniki obliczeń, tzn. średnią z 10 minimalny h długości trasy komiwojażera uzy- skanych po wykonaniu n iteracji. Rozważamy różne warianty parametru n odpowiadające ostatniej kolumnie tabeli ze str. 123.

Ostatnia kolumna tabeli zawiera dla porównania długości trasy komiwo- jażera uzyskane algorytmem Volgenanta i van den Ilouta (VOL).

m n MC SA: 20PT SA:RANPER VOL

15 230 5.087 3.877 5.050 3.281

299 5.013 3.716 4.978

2302 4.594 3.316 4.674

2995 4,588 3.337 4.730

20 230 7.395 5.258 7.407 3.984

299 7.351 4.960 7.369

2302 6.683 4.108 6.731

2995 6.750 4.106 6.773

m n MC SA: 20PT SA: RANPER VOL

50 230 22.596 15.037 22.302 5.622

299 22.208 13.715 21.963

2302 21.305 6.872 21.274

2995 21.361 6.357 21.226

Z powyższych obliczeń widzimy, że losowanie typu 20PT ujawnia swą przewagę dla algorytmu SA, ponadto przyspiesza obliczenia. Jednak istotne różnice powstają dopiero dla bardzo dużej liczby iteracji. W przeciwieństwie do rezultatów uzyskanych dla funkcji testowych f

i }

rozdziale 2.5 za- uważamy, że metoda MC zachowuje się o wiele gorzej w zadaniu komiwoja- żera. Wyniki dla algorytmu MC są zbliżone do wyników algorytmu SA z lo- sowaniem sąsiedztwa typu RANPER (ale w badanych przypadkach zawsze gorsze), co wskazuje na duże znaczenie odpowiednio wybranej struktury otoczeń w problemach kombinatorycznych. Przeprowadzone eksperymenty sugerują również postawienie problemu charakteryzacji zadań, w których algorytm MC byłby lepszy od algorytmu SA.

3. Przypadek z zaburzeniami

W dalszej części pracy rozważać będziemy sytuację, w której badana funkcja poddana jest pewnym losowym zaburzeniom. W praktycznych za- daniach optymalizacji dyskretnej często pojawia się problem zaburzeń war- tości funkcji. Wynika to na przykład z faktu niedokładności w pomiarach badanej wielkości lub też ze zbyt dużego kosztu dokładniejszych pomiarów.

W teoretycznych badaniach algorytmów optymalizacji musimy uwzględnić wpływ zaburzeń na takie własności stosowanych algorytmów jak zbieżność i efektywność osiągania ekstremum globalnego.

3.1. Prosta metoda Monte Carlo

Wyprowadzimy teraz pewne oszacowanie dla miary precyzji algorytmu pro- stej metody Monte Carlo, tzn. dla wielkości P(An(f*)) w przypadku, gdy występują losowe zaburzenia wartości funkcji.

Niech Fe oznacza dystrybuantę zaburzenia , natomiast G dystrybuantę zmiennej losowej (patrz rozdz. 1.4, str. 5).

Przyjmijmy również następujące oznaczenia:

fi* = mm \fk - fi\, S* = max \fk - fi|, a = ^{ć *} o

Na początku pracy założyliśmy różnowartościowość funkcji /, skąd mamy

fi* > 0.

Algorytmy stochastyczne w optymalizacji dyskretnej

Parametry a i fi dla danej funkcji / i ustalonego parametru skali cr rozkładu zaburzeń stanowią, wskaźniki względnego poziomu zaburzeń.

L

3. Zachodzi następująca nierówność i - 1

p(Avif,: » > p ‘ («) wH{n, N) gdzie d — [Ni ]\,

+oo

(//. N) dla N — d — 1 > /?, dla N — d — 1 / n

w

(n, N) = j [1 — G(t — er)]**•[ 1 — G{t - [n - \_)a)]dG{t).

D ow ód. Bez zmniejszania ogólności założymy, że jj < h < • • • < I

oraz L < /2 < . . . < in.

Zdarzenie zachodzi na pewno wtedy, gdy zajdą jednocześnie dwa następujące zdarzenia:

1 ) wylosowany podzbiór {/ 1 , ..., in] zawiera przynajmniej jeden element ze zbioru { I,. . ., d -f 1 } oraz

2 ) element o najmniejszej zaburzonej wartości funkcji spośród elemen- tów

,..., im oznaczany przez jest dokładnie tym samym ele- mentem, który zostałby wyznaczony przy obserwacjach funkcji bez zaburzeń ( czyli ?'* = ii).

Oba te zdarzenia są, na mocy założeń na początku pracy, niezależne. Zatem prawdopodobieństwo zdarzenia /l7?(/*) jest nie mniejsze od iloczynu praw- dopodobieństw 7 tych dwóch zdarzeń. Prawdopodobieństwo pierwszego z nich wrynosi, jak poprzednio (przypadek bez zaburzeń, w7zór ( 1 ))

])(n, iV, d) r t 1)

o ’ dla N — d — 1 > n l 1 dla N — d — 1 < ii a prawdopodobieństwu) drugiego oznaczymy przez w(n,N).

Zatem zachodzi nierówność

P ( M K ) ) >

N — d —l

W{ 7 7., N )

) )w(n, N) dla, N — d — 1 > z?, dla N — d - 1 < n.

Pozostaje jeszcze obliczenie w{n,N). Na mocy powyższych rozważań mamy

M n , N ) = p ( hi = h ) = ^ ( m i n { fh , • • ■, f i n } = h 1 )•

Zauważmy, że zachodzi oszacowanie P(i*n = ń ) = P( fh < f ik,k =

^ P[£il ^ + ń*

2 ,..., 77.) = P{{eik - eh ) > - f ik - fi,))

€i 1 < e«3 + 2<5*, . . . , < €jn + (77. — !)<*>*)

= J ^dFe {xi)---(lFe {xn)

{x1<x2 + &*,-..,Xi<xn+(n-l)5*}

= f [1

- F ,{x - 6,)] • • • [1 - F t(x - (n - 1 )6 , )}dF (x)

— oo

-\-oo

- - j [1 — G ( i — a)] • • • [1 — G it — (n — 1 )a )}r lG ( t) = wH(n, N) .

— OO

S tąd uzysku jemy w(n,N) > wH(n,N) i tezę tw ierdzen ia . ■ D la dużych parame trów a i (3 , czy l i d la „ma łych” zaburzeń , oszacowa- n ie z poprzedn iego lema tu da je wyn ik co raz b l iższy w ie lkośc i P(Av(f*)) z przypadku bez zaburzeń . Z jaw isko to , k tóre można ok reś l ić jako c iąg łość ze wzg lędu na zaburzen ia , wy jaśn ia nas tępu jący lema t .

L

4 . Zachodz i

l im wH(n , N) = 1 .

a —roo

Dowód wyn ika z de f in ic j i w ie lkośc i wys tępu jących w lemac ie oraz z fak tu , że

l im G ' ( : r} = 0 ,

x —► — OO

k tóry wyn ika z ko le i z w łasnośc i G jako dys trybuan ty rozk ładu prawdopo-

dob ieńs twa . ■

Podobn ie jak w sy tuac j i bez zaburzeń chcemy dobrać taką war tość ne ( N , i] , 7 ) , aby zapewn ić odpow iedn ią jakość rozw iązan ia tzn . aby d la da- nych r j oraz 7 spe łn iona by ła n ierówność

p{^rA fn)) > 7 gdy n > ne (N, ? / , 7 ) .

Na mocy lema tu 3 n ,(N , / / , 7 ) jes t (na jmn ie jszym) rozw iązan iem n ierów- nośc i

(7 ) p { ? ? - , A r , d)wH (n, N) > 7 .

Na pods taw ie powyższe j n ierównośc i n ie można podać jednoznaczne j regu ły pozwa la jące j na wyznaczan ie w ie lkośc i ne (N , r / , 7 ) . W w ie lu kon- kre tnych przypadkach n ierówność (7) n ie ma rozw iązan ia , gdyż w ie lkość wH(n,N) może być d la wszys tk ich war tośc i n mn ie jsza od 7 . W te j sy- tuac j i nasuwa , s ię nas tępu jący pomys ł mody f ikac j i a lgory tmu m in ima l izac j i po lega jący na . tym , aby w każdym kroku a lgory tmu , czy l i w każdym punk- c ie i G 1 uzysk iwać n ieza leżn ie k war tośc i funkc j i {// . ,• •• , / /) .} (schema t w ie lokro tnych pom iarów) ; do da lszych ob l iczeń pods taw iamy uśredn ione war tośc i

k

Algorytmy stochastyczne w optymalizacji dyskretnej

Procedura uśredniania pozwala zmniejszyć wpływ zaburzeń, dzięki temu.

ze

Var + ••• + * -Varc.^ < Varcjj Stosując nierówność Czebyszewa marny

V i>0 P + ... + c

< t \ > 1 - fc?Vare"

Zatem dla ustalonego poziomu zaburzeń można dobrać współczynnik powtórzeń k tak, aby wielkość tvH[n,N) była bliska jedynce.

Przeprowadzimy teraz analizę nierówności (7) dla dużych N . Z poprzed- niego rozdziału dotyczącego przypadku bez zaburzeń wiemy, że

lim p{n, N, d) = 1 — (1 — ?/)n.

jV-+oo

Powyższe rozważania prowadzą do następującego schematu znajdowania rozwiązania nierówności (7):

1) Dobieramy współczynnik krotności obserwacji k tak, aby wH[n, N ) > 7 .

Wielkość wH(n, iV) spełniającą powyższą nierówność oznaczmy przez

■w 0.

2 ) Jako rozwiązanie przyjmujemy

ne( N, ?/, 7 ) = r v ln(l - 3-)

■ ln( 1 — ?/) i-

Przyjmując konkretne rozkłady zaburzeń uzyskamy konkretne wartości liczbowe. Rozważmy dwa przypadki, najczęściej stosowane do opisywania zaburzeń występujących w realnych zadaniach.

1. Zaburzenia < 7 . mają rozkład normalny tzn. Fa ~ jV(0,<r1 2). Wtedy uzyskujemy wzór

+ CO

wH(n,N) = j 4>{t)[l - J>( / - o)] ■■•[!- J>( / - ( 77 . - 1 .)cxj]dt, gdzie $ oznacza dystrybuantę, a (f> gęstość standardowego rozkładu normalnego iV( 0 , 1 ).

II. Zaburzenia ep mają rozkład jednostajny tzn. Fa ~ U {—er, a). Wtedy odpowiedni wzór ma postać

1 J-

wH{ ¹¹ , N) = - j [1 — G(t — «)] •••[! — G{t — [n — 1 )o )]<■//.

- 1

gdzie G jest dystrybuantą rozkładu jednostajnego U{ — 1 , 1 ).

W poniższej tabeli przedstawiamy przykładowe wartości 7 ) dla dwóch rodzajów zaburzeń.

Zaburzenie a 7 fi k ne{N, i], 7 ) normalne 1.0 0.9 0.01 4 380

0.5 0.9 0.01 14 512 jednostajne 1.0 0.9 0.01 2 397 0.5 0.9 0.01 5 422

Z podanych wielkości odczytujemy np., że dla zaburzenia normalnego z a = 1 należy wykonać 4-380 = 1520 obliczeń funkcji, aby uzyskać gwarancję podobnej jakości rozwiązania jak w przypadku bez zaburzeń (w którym wystarczy tylko 230 obliczeń funkcji).

3.2. Algorytmy sekwencyjne

Zmodyfikujemy schemat (A) algorytmów sekwencyjnych z rozdziału 2.2, uwzględniając występowanie losowych zaburzeń. Na każdym etapie algo- rytmu zamiast wartości dokładnych obserwujemy

fi — fi 4~ €ii

gdzie Ci jest zmienną losową odpowiadającą zaburzeniu. Podamy teraz kon- strukcję ciągów zmiennych losowych (Yk,k > 1 ). Zmienna Yk modeluje stan optymalny na k -tym etapie algorytmu, natomiast zmienna llfi. jest nie- obserwowalną zmienną losową, która modeluje zaburzenie losowe, z jakim obliczana jest na k-tym etapie algorytmu różnica wartości funkcji. Zmienna

\\'\ ma postać:

Wk = €j - Ci.

Niech Yi będzie wybrane w sposób dowolny. Niech Fk oznacza dalej dys- trybaantę zmiennej losowej Wk z ciągu (Wk,k > 1). Dla danych elemen- tów !'i, ..., Yk„ następująco definiujemy zmienną losową Yk+ \: według pew- nego ustalonego rozkładu prawdopodobieństwa {gyk,j, j G A'} losujemy stan j G A' i wtedy

f i gc|y fj - frk > 0 i Bk < oi(fj ~ fyk, Tk) [B] Yk+1 = < j gdy fj - f Yk < 0 i R k < [3(fj - fyk,Tk)

1 Yk w przeciwnym przypadku Obliczymy następujące prawdopodobieństwo

(8) Pr{Yk+1 = j |1 j , .. ., n _ ! , n = i, Wu . . . y W k-i, Wk I (JijBuWi fj - fi + A, Tk)}, jeśli fj - fi Y A > 0 1 (JijFu{ft{f3 - fi + A, Tk)}, jeśli fj - fi Y A < 0

= A }

V?;, j G A, A G R

gdzie Fu jest dystrybuantą rozkładu jednostajnego na [ 0 , 1 ].

Algorytmy stochastyczne w optymalizacji dyskretnej

Z konstrukcji zmiennych losowych Y\. i wzoru 8 wynika, że ciąg (Yk. k >

1 ) jest łańcuchem Markowa z przestrzenią stanów X i macierzą przejść w jed- nym kroku P{k) = [P.i:j{ k)]ijjeX, gdzie

Pij{k) = Pr{I'Ar+i = j\Yk = i] = Ewk {Pr{l*+i = j\Yk = i, Wk}}

= J [J i j Eu { ot (fj - fi + A, Tk)} dFk (A) {\:\> fi — f j }

+ f a . j f ' u W f j - f. + K T k)}dFk{\)

= gij(iij{Tk),

« ( T k ^{) =} J Fu { ot {fj - f t + A, Tp)} d Fj-(A)

{A: A

+ / F u W f j - f i+ X,Tk)}dFk(\)

{A: A < / , - / , }

Algorytmy typu (B) stanowią uogólnienie schematu (A) na przypadek funkcji obserwowanych z zaburzeniami losowymi.

W dalszym ciągu rozważań uprościmy ogólny schemat (B) wprowadzając następujące założenie:

Zmienne ek zaburzające wartości funkcji są niezależnymi zmiennymi lo- sowymi o nośniku ograniczonym z dołu. Wtedy również zmienne Wk są ograniczone z dołu.

3.3. Zbieżność ciągu rekordów algorytmu

Oprócz ciągu {Yk,k > 1) interesować nas również bedzie ciąg tzw. rekordów algorytmu (B), określony w sposób następujący:

rk = arg min f(Ym)

1 <.m<.k

Gdy występują zaburzenia losowe, to ciąg (rk,k > 1 ) jest ciągiem rekor- dów estymowanych rekurencyjnie na podstawie przybliżonych informacji o wartościach minimalizowanej funkcji.

Udowodnimy następujące twierdzenie.

Tw i e r d z e n i e 4.

Jeśli istnieje liczba naturalna ko taka, że ciąg {a(k), k >

ko) jest nierosnący, to z rozbieżności szeregu (LM{kM) wynika zbież- ność ciągu rekordów (rk,k > 1) 2 : prawdopodobieństwem 1 do zbioru X*.

Dowód. Niech d = min.{#.;j : gtl > 0}.

Wtedy mamy

Pij{k) > d a{k), Vk>ko V{i,jeX:gij>o}