inynierskie metody numeryczne D 10230 bszafranagh edu pl

inżynierskie metody numeryczne D 10/230, bszafran@agh. edu. pl konsultacje: środy 8: 30 -10: 00 http: //galaxy. uci. agh. edu. pl/~bszafran/imn 10 cel przedmiotu: przygotowanie do pracy w zakresie numerycznego modelowania zjawisk i urządzeń stosowanego w zagadnieniach techniki (inżynierii) i nauki W technice: inżynieria obliczeniowa: modelowanie i symulacja zjawisk i urządzeń. badania i optymalizacji procesów produkcyjnych oraz produktów. Obliczenia naukowe – interpretacja i przewidywanie danych doświadczalnych, zrozumienie obserwacji, przewidywanie nowych zjawisk. Modelowanie naukowe/inżynieryjne: metody podobne, różnica w celu oraz obiekcie badań

inżynierskie metody numeryczne Postęp w technologii komputerowej i metodach obliczeniowych sprawia, że znaczenie inżynierii obliczeniowej w przemyśle i nauce ciągle rośnie - co dotyczy wszystkich gałęzi przemysłu i dziedzin naukowych. Personel wykwalifikowany do prowadzenia obliczeń jest na świecie coraz bardziej poszukiwany. Najczęściej poszukiwani są specjaliści z danej dziedziny inżynierii obeznani z podstawami modelowania komputerowego - stosujący gotowe pakiety obliczeniowe. Wykład dla kierunku informatyka stosowana – nie pracujemy z pakietami obliczeniowymi, metody implementujemy od podstaw – przygotowanie do działalności w opracowaniu oprogramowania obliczeniowego.

Metody wyznaczania własności układów: 1) Metody teoretyczne (ścisłe, analityczne) ograniczone do prostych problemów (w nauce – te akurat są często najważniejsze) lub wyidealizowanych modeli. Idealizacja oparta na intuicji, które bywa błędna. 2) Badania doświadczalne nieodzowne i najważniejsze, ale drogie (i / lub czasochłonne) często same w sobie nie pozwalają na zrozumienie zjawiska przydatne wsparcie ze strony obliczeń ścisłych lub przybliżonych 3) Symulacje numeryczne pozwalają na rozwiązywanie dokładnych równań z kontrolowalną dokładnością często do wykonania taniej i szybciej niż badania doświadczalne pozwalają prześledzić wyniki w funkcji dowolnych parametrów – pełna informacja o możliwych do osiągnięcia własnościach

Tematyka wykładu: głównie rozwiązywanie równań różniczkowych: zwyczajnych i cząstkowych. Równania różniczkowe – opis zjawisk wprowadzony w XVIII / XIX w. Możliwy w praktycznym zastosowaniu w modelowaniu od niedawna – od pojawienia się wydajnych komputerów [np. Równania dynamiki płynów znane od połowy XIX wieku Navier/Stokes - dopiero od niedawna (ok. 30 lat) są rozwiązywane poza najprostszymi przypadkami] Metody numeryczne – przybliżone i wydajne rozwiązania równań równie stare jak same równania, ale przed komputerami stosowalne tylko do najprostszych problemów (co nie znaczy, że bez znaczenia – odkrycie Neptuna)

Znaczenie modelowania numerycznego rosło i będzie rosło z rozwojem sprzętu obliczeniowego Rok 1949 EDSAC (lampowy) 1997 ASCI Red (symulator eksplozji jądrowych) 2002 NEC Earth Simulator (modelowanie klimatu) 2009 IBM Blue Gene / Q FLOPS 102 1013 1015 Pamięć 2 k. B 300 GB 10 TB 500 TB Metody obliczeniowe również się rozwijają (za M. Schaeferem, computational engineering) Rok 1970 1975 1980 1985 1990 2005 Eliminacja Gaussa Metoda Gaussa-Seidla Nadrelaksacja Metoda gradientów sprzężonych Metody wielosiatkowe Siatka adaptowana Tempo rachunków 1 15 250 1 k 5 k 50 k

Symulacje numeryczne są ze swej natury interdyscyplinarne (matematyka, metody numeryczne, nauki ścisłe, konkretna dziedzina inżynierii / nauki + programowanie) - gdzie trzeba będę starał się podawać elementarną wiedzę z zakresu fizyki opisywanych zjawisk. Matematyka numeryczna Informatyka Symulacja numeryczne Fizyka/Chemia Dziedzina pochodzenia problemu (inżynieria/nauka) Tylko w ujęciu inter – symulacje są użyteczne (interesujące)

Miejsce numeryki w rozwiązywaniu problemów Problem (naukowy/inżynieryjny) Rozwiązanie Dane doświadczalne, modele matematyczne Analiza i interpretacja Równania różniczkowe / warunki brzegowe Weryfikacja i korekta modelu Generacja siatki, dyskretyzacja (czasu / obszaru całkowania) * Układy równań algebraicznych wartości użyteczne / mierzalne Przetworzona informacja *** Obróbka danych *** Algebraiczne algorytmy numeryczne ** programowanie *** *wykład (FDM, FVM, FEM, BEM) ** wykład (ten lub KSN) *** laboratorium Rozwiązanie numeryczne (milion liczb) ***

Szczegółowy program wykładu: galaxy. uci. agh. edu. pl/~bszafran/imn 10/planw 10. pdf Literatura: Press, Numerical Recipes (The art of scientific computing). Haupt, Practical Genetic Algorithms. Weinberger, A first course in partial differential equations. Koonin, Computational Physics. Solin, Partial Differential Equations and the Finite Element Method. Zienkiewicz, Finite Element Method, its basis & fundamentals. Lienhard, A Heat Transfer Textbook. Sabersky, Fluid flow : a first course in fluid mechanics. Quarteroni, Numerical mathematics. Trefethen, Finite difference and spectral methods for ordinary and partial differential equations. Cichoń, Metody Obliczeniowe. Sewell, The numerical solution of ordinary and partial differential equations. Evans, Numerical Methods for Partial Differential Equations R. Grzymkowski, A. Kapusta, I. Nowak, D. Słota, Metody Numeryczne, Zagadnienia Brzegowe Schafer, Computational Engineering, Introduction to Numerical Methods

Laboratorium: staramy się aby związek wykładu z laboratorium był bliski 1: 1 Tematy (za treść odpowiedzialny – wykładowca) - na stronie http: //galaxy. uci. agh. edu. pl/~bszafran/imn 10 Podane jest w tej chwili około 10 tematów – ich treść może się zmieniać, ale nie później niż 7 dni przed zajęciami laboratoryjnymi. 4 Grupy: Wtorek 14: 30 -17: 30 Środa 18 -21 Środa 14: 30 -17: 30 Środa 18 -21 dr inż. Małgorzata Krawczyk dr inż. Przemysław Gawroński dr inż. Krzysztof Malarz dr inż. Maciej Wołoszyn

Ocena z laboratorium: Średnia arytmetyczna z aktywności i raportów. Aktywność: oceniana na podstawie wyników i źródeł przesłanych prowadzącemu pod koniec zajęć. Raporty: pełne, lecz możliwe krótkie opisy uzyskanych wyników. Aby raport kwalifikował się do zaliczenia (dt): musi zawierać poprawne wyniki Raport dobry – ma wyniki przedstawione w sposób czytelny i jasno opisany Bardzo dobry raport musi zawierać samodzielnie napisany zwięzły wstęp i wnioski. Raport należy przysłać prowadzącemu w terminie 14 dni (w semestrze zimowym 7 dni) od zakończenia zajęć. Raporty przysłane później nie będą czytane/oceniane.

Ocena z laboratorium: Student ma prawo do dwukrotnego poprawiania zajęć laboratoryjnych. Dwa dodatkowe zajęcia laboratoryjne. Na każdym z nich można poprawić jedną aktywność i jeden raport albo anulujemy każdemu studentowi dwie najgorsze oceny z aktywności / raportu

Egzamin: pisemny (test), 5 zadań do rozwiązania. Przystąpić mogą osoby z zaliczonymi obydwoma semestrami lab. Ocena: 1/3 średniej z zaliczeń + 2/3 wyniku testu Dla potrzeb egzaminu średnia z zaliczeń policzona będzie z górnych widełek dla danej oceny: Regulamin AGH: 91 – 100% bardzo dobry (5. 0); 81 – 90% plus dobry (4. 5); 71 – 80% dobry (4. 0); 61 – 70% plus dostateczny (3. 5); 50 – 60% dostateczny (3. 0); (3. 5 + 5)/2 = (70%+100%)/2 = 85% Przykład: osoba z dwoma czwórkami rozwiązuje 2 zadania na 5 1/3 80% + 2/3 40% >50 % Wniosek: dobre zaliczenie bardzo pomaga zdać egzamin

Zwolnienie z egzaminu: Student, którego średnia z zaliczeń laboratoryjnych będzie odpowiadała 4. 5 (>80%) - może zostać zwolniony z egzaminu [średnia z zaliczeń zostaje wpisana wtedy jako wynik egzaminu]: Pary ocen uprawniających do zwolnienia z egzaminu: 4 + 4. 5 4 + 5 3. 5 + 5 4. 5 + 5 → 4. 5 → 5

Obecny plan laboratorium: Zajęcia 1: algorytm kolonii mrówek Zajęcia 2: symulowane wygrzewanie ( Środki przeszkód do zadania 2). Zajęcia 3: wariacyjne MC z algorytmem Metropolisa dla równania Poissona Zajęcia 4: schematy jawne i niejawne dla rrz Zajęcia 5: ekstrapolacja Richardsona: szacowanie błędów i automatyczny dobór kroku czasowego Zajęcia 6: problemy sztywne Zajęcia 7: niejawne metody RK Zajęcia 8: problem brzegowy 1 D, metoda Numerowa Zajęcia 9: iteracja wielosiatkowa dla równania Laplace'a Zajęcia 10: metody relaksacyjne Zajęcia 11: przepływ potencjalny Zanim dojdziemy do równań różniczkowych – mniej twarde rachunki, być może mniej ważne – ale nie mniej ciekawe niż główny wątek wykładu. . . Zaczynamy od mrówek. . .

optymalizacja kombinatoryczna – dla zmiennej dyskretnej Przykład: najkrótsza sieć (kanalizacyjna, energetyczna, światłowodowa) dla miast powiatowych województwa małopolskiego:

optymalizacja kombinatoryczna – dla zmiennej dyskretnej Przykład: najkrótsza sieć (kanalizacyjna, energetyczna, światłowodowa) dla miast powiatowych województwa małopolskiego: Województwo zakodowane w postaci grafu – poszukiwane najkrótsze drzewo spinające. 31

optymalizacja kombinatoryczna – dla zmiennej dyskretnej Przykład: najkrótsza sieć (kanalizacyjna, energetyczna, światłowodowa) dla miast powiatowych województwa małopolskiego: Województwo zakodowane w postaci grafu – poszukiwane najkrótsze drzewo spinające. 31 Rozwiązanie:

optymalizacja kombinatoryczna – dla zmiennej dyskretnej Przykład: najkrótsza sieć (kanalizacyjna, energetyczna, światłowodowa) dla miast powiatowych województwa małopolskiego: Województwo zakodowane w postaci grafu – poszukiwane najkrótsze drzewo spinające. 31 Problem równie łatwy jak regresja liniowa Rozwiązanie: rozwiązanie dane przez algorytm zachłanny Kruskula : tworzymy las dodając po kolei najkrótsze krawędzie tak aby nie utworzyć pętli: dostaniemy najlepsze rozwiązanie a. zachłanny: kieruje się regułą maksymalnego zysku w najbliższym kroku (w szachach sukcesów nie odnosi)

Limanowa-Nowy Sącz C 41 T 43 42 B 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS K NT

Wadowice-Chrzanów C 41 T 43 42 B 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS K NT

Myślenice-Kraków C 41 T 43 42 B 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS K NT

Bochnia-Limanowa C 41 T 43 B 42 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS K NT

Wadowice-Myślenice C 41 T 43 42 B 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS K NT

Nowy Sącz-Gorlice C 41 T 43 42 B 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS K NT

bezpośrednie połączenie CK już się nie przyda Kraków-Bochnia C 41 T K 42 43 B 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS NT

M-NT B-T C 41 T K 42 43 B 31 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS NT

31 Algorytm zachłanny – skuteczny więc - problem najkrótszego drzewa spinającego jest łatwy. Rozwiązanie: Złożoność dla najlepszej implementacji O(|V|log|V|), V – liczba wierzchołków

Problem najkrótszej drogi z wierzchołka A do wszystkich pozostałych w grafie (przeszukiwanie grafu wszerz z oznaczaniem wierzchołków) Każda krawędź ma wagę (długość) – ustalone i definiujące problem Wierzchołkom przyznajemy: Wagi: najkrótsza znaleziona odległość z wierzchołka A Koloru: biały oznacza, że waga tymczasowa, czarny, że ustalona (krócej nie będzie) Etykiety: którędy do A Wagi, kolory i etykiety wierzchołków zmieniają się w czasie działania algorytmu 0) Oznacz wszystkie wierzchołki kolorem białym. Przypisz wierzchołkowi startowemu (A) wagę 0 1) 2) Znajdź i zaczerń biały wierzchołek v o najmniejszej wadze Dla każdego B - białego wierzchołka przyległego do v: jeśli waga B > waga v + waga krawędzi (v, B) lub jeśli waga B nie została wcześniej zdefiniowana waga B: =waga v + waga krawędzi (v, B) oraz ustaw etykietę B na v Jeśli są jeszcze białe wierzchołki idź do 1 3) (złożoność V 2)

przykład: najkrótsze trasy z Gorlic do pozostałych miast 1) Gorlice malujemy na czarno, miastom sąsiednim nadajemy wagi – odległości od Gorlic i indeks G. C 41 T 45 G 42 B 44 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS 36 G NT K 3) sąsiedzi Nowego Sącza otrzymują próbne wagi i etykiety C 41 T 45 G 42 B 44 31 36 M 42 45 42 32 W L 59 NS 65 43 23 S 36 G 57 NS 36 G NT 93 NS K 2) Szukamy białego miasta o najmniejszej wadze i malujemy je na czarno (Nowy Sącz), wagę czarnego miasta ustalamy (mniejszej nie będzie) C 41 T 45 G 42 B 44 31 36 M 42 45 42 32 W L 65 43 23 S 36 G 57 NS 36 G NT K

4) Najmniejszą wagę ma teraz Tarnów, C C 41 5) Następnie Limanowa 41 89 T T 45 G 42 B 44 31 36 M 42 45 42 32 W L 59 NS 65 43 23 S 36 G 57 NS 36 G NT 93 NS K 6) Bochnia C 41 131 B 89 T T 45 G K 42 B 44 31 36 M 42 131 B? 45 42 32 W L 59 NS 101 L 65 43 23 S 36 G 57 NS 36 G NT 93 NS do Myślenic jednak bliżej przez Limanową 89 T T 45 G 42 B 44 31 do Bochni z Gorlic 91 L ? 36 M 42 45 bliżej przez Tarnów niż 42 32 W L 59 NS przez Limanową 101 L 65 43 23 S 36 G 57 NS 36 G NT 93 NS K 7) Po Nowy Targu – Myślenice, z nich bliżej do Wadowic

ostatecznie Np. : z Chrzanowa do Gorlic trafimy po etykietach Zamiast stosować algorytmu można zrobić model z nitek i koralików, potem naciągnąć koraliki oznaczające Chrzanów i Gorlice

Widzieliśmy, że dwa ważne problemy mają efektywne, deterministyczne, dokładne rozwiązanie Niektóre problemy są jednak obiektywnie trudne (nie istnieje algorytm o złożoności wielomianowej): wybór najkrótszej zamkniętej trasy przez wszystkie miasta (problem komiwojażera): Odwiedzić wszystkie miasta w cyklu zamkniętym w takiej kolejności aby pokonana trasa była najkrótsza. -algorytm deterministyczny rozwiązujący problem dokładny z wielomianową złożonością nie istnieje, gdy problem o dużym rozmiarze należy rozwiązać – stosuje się podejścia heurystyczne lub Monte Carlo Algorytm zachłanny dla komiwojażera: Klasyczny problem testowy dla algorytmów optymalizacyjnych ruszaj do najbliższego miasta, którego jeszcze nie odwiedziłeś. - rozsądny: wyeliminuje przynajmniej długie przejazdy bez zatrzymywania się

Zachłanne rozwiązanie jest optymalne (choć nie najgorsze) Rozwiązanie zachłanne: start ze Szczecina: PL: 46 miast Najlepsze Szukana jest permutacja - przejrzeć wszystkie N! - niewykonalne 46!=55026221598120889498503054288002548929616517529600000 najlepszy algorytm dokładny: złożoność rzędu 2 N –lepiej niż n!, ale wciąż zbyt wiele 246=70368744177664 Gdy problem zbyt trudny by go rozwiązać dokładnie przy pomocy algorytmu deterministycznego – można zadowolić się przybliżonym (heurystycznym) lub próbować je poprawić przy pomocy MC

Problem obiektywnie trudny = gdy najlepszy deterministyczny algorytm nie zakończy swojego działania w skończonym czasie klasy złożoności obliczeniowej Problemy decyzyjne: z odpowiedzią tak/nie Problemy NP P NP-zupełne Schemat obowiązuje pod warunkiem że P NP NP – można sprawdzić odpowiedź w czasie wielomianowym zadanie rozkładu na czynniki pierwsze liczby 136117223861 nieznany jest wielomianowy algorytm (na komputer klasyczny) ale jeśli ktoś nam poda odpowiedź 104729 1299709 - szybko sprawdzimy. P – problemy, w których istnieje algorytm o wielomianowej złożoności (nie ma dowodu, że P NP. ) NP – zupełne (najtrudniejsze) – można do nich sprowadzić dowolny problem z NP z nadkładem wielomianowym. Jeśli jeden z problemów NP. -zupełnych zostanie rozwiązany w czasie wielomianowym, to P=NP.

Problemy NP P NP-zupełne F Faktoryzacja jest na pewno NP, wydaje się, że nie jest P i że nie jest NP-zupełna. [„Wydaje się, że nie P” na tyle bardzo, że protokół klucza publicznego RSA oparty na niewykonalności faktoryzacji w krótkim czasie] NP. -zupełne: problem spełnialności binarnego układu logicznego, problem komiwojażera, izomorfizmu grafów, kliki, kolorowania wierzchołków grafu i inne. W praktyce problemy, które nie są P – stają się niemożliwe do dokładnego rozwiązania dla dużych rozmiarów zadania

najkrótsza trasa z A do B – problem łatwy (bo wielomianowy algorytm znany) najkrótsza trasa po wszystkich miastach – problem trudny (bo algorytm wielomianowy nieznany) Inna znana para pozornie podobnych problemów o skrajnie różnej złożoności obliczeniowej: problem istnienia cyklu Eulera i cyklu Hamiltona Cykl (zamknięta ścieżka) Eulera Zadanie: zaplanować trasę spaceru: przejść po każdym moście dokładnie raz i wrócić do punktu wyjścia (odwiedź wszystkie krawędzie grafu dokładnie raz)

Cykl Eulera w grafie istnieje wtedy i tylko wtedy gdy wszystkie jego wierzchołki są stopnia parzystego 3 stopień wierzchołka = liczba przyległych krawędzi 3 5 3 Jeśli istnieje to parzyste stopnie Jeśli parzyste stopnie to można wskazać cykl E przy każdym przejściu przez wierzchołek używamy 2 krawędzi zaczynamy spacer od dowolnego wierzchołka usuwając z grafu przebyte krawędzie, wrócimy do wierzchołka startowego bez rozspójniania grafu

Cykl Hamiltona odwiedzić wszystkie wierzchołki bez powtarzania przemarszu przez żadną z krawędzi wracając do punktu wyjścia (nie trzeba korzystać z każdej krawędzi, co wcale nie ułatwia zadania, bo problem istnienia CH – np. -zupełny ) graf planarny (rzut środkowy dwunastościanu)

cykl Hamiltona dla dwunastościanu

Jeśli wiemy, że problem NP-zupełny, a rozmiar problemu duży – poszukajmy rozwiązania przybliżonego Metoda dokładna nie zadziała w skończonym czasie. Jeśli nie wiemy jak szukać- poszukajmy losowo. Lecz: Całkiem ślepe przeszukiwanie losowe nie różni się od przeglądania wszystkich rozwiązań: prawdopodobieństwo znalezienia najlepszego jest żadne, a i rozsądnego znikome. Problem komiwojażera dla 20 miast w pd-wsch Polsce Wszystkich permutacji jest 20!=2432902008176640000. Najlepsza trasa znaleziona po 1000 prób (długość 89. 12 [j. umowne] ) Widać, że kiepska: 1) skrzyżowane trasy 2) krócej będzie Tarnów-Nowy Sącz-Kraków Katowice

Najlepsza trasa znaleziona po 1000 losowaniach (długość 89. 12 [j. umowne] ) Algorytm zachłanny start z Częstochowy 68. 73 Wniosek: jeśli przeszukiwanie MC ma być skuteczne – nie może być całkiem losowe

Optymalizacja wg podejścia „kolonii mrówek” Dorigo M, Gambardella LM BIOSYSTEMS 43 73 1997 Mrówki zostawiają na swojej trasie szlak zapachowy (feromonowy). Pozostałe mrówki starają się go śledzić. Ręcznie położony szlak zapachowy Ślad zapachowy: 1) wyznacza trasę między mrowiskiem a źródłem pożywienia 2) pozwala na jej optymalne wytyczenie 3) umożliwia adaptację do zmieniającego się środowiska źródło rysunków: Eberhart, Shi, Kennedy: Swarm Intelligence

źródło rysunków: Eberhart, Shi, Kennedy: Swarm Intelligence Dodatnie sprzężenie zwrotne (najsilniejszy ślad na najkrótszej drodze).

Szlak zapachowy i adaptacja trasy dla zmienionego otoczenia wstawiona przeszkoda losowo na prawo/lewo ? ? Na krótszej drodze szlak feromonowy szybciej domknięty Dodatnie sprzężenie zwrotne Dorigo M, Gambardella LM BIOSYSTEMS 43 73 1997

Algorytm kolonii mrówek dla problemu komiwojażera Założenia: 1 2 3 4 5 W metodzie: grupa przemieszczających się wędrowców symulujących zachowanie mrówek Preferowany szlak o największym nasileniu feromonu Wzrost feromonu na najkrótszej trasie Komunikacja między mrówkami przez szlak feromonowy Strategia najbliższego sąsiada włączona w algorytm (mrówki kierują się nie tylko feromonem, ale również starają się dostać do miast raczej bliżej niż dalej położonych) Sztuczne mrówki potrafią więcej niż prawdziwe: (1) Pamiętają swoją trasę (2) Ustalają, która przeszła najkrótszą drogę

Miasta połączone w graf pełny: krawędź między wierzchołkami i oraz j opisana (1) długością d(i, j) stała, definiuje problem (2) natężeniem feromonu f(i, j) zmienia się w trakcie nauki Z długości (stałej) i feromonu (zmiennego) budujemy funkcję preferowanego ruchu: t(i, j) = f(i, j) / [d(i, j)]b im większe t(i, j) tym większe będzie pstwo przejścia z i do j b– parametr [znaczenie feromonu/długości krawędzi] b=0 tylko feromon, duże b -tylko długość optymalne: b=2

Mrówki wykonują wiele obiegów przed każdym z nich są losowo rozmieszczane po miastach przykład:

Wybór mrówki: jest w mieście i do którego miasta iść? Każda krawędź opisana przez t(i, j) ? i Tu nie pójdzie bo już była mrówka kieruje się doświadczeniem zdobytym przez mrowisko, nawet gdy wybór miasta następuje przy użyciu generatora liczb losowych Parametr q 0 z przedziału (0, 1). q 0 – jeśli wysokie mrówki bardziej skłonne korzystać z doświadczenia społeczności – jeśli niskie mrówki chętniej eksperymentują z nowymi trasami Losujemy q z przedziału (0, 1): jeśli q<q 0 [doświadczenie] mrówka idzie do nie odwiedzonego jeszcze miasta j, dla którego t(i, j) = max jeśli q>q 0 [eksploracja] optymalne q 0=0. 9 miasto j jest losowane z rozkładem, w którym prawdopodobieństwo jest proporcjonalne do t(i, j)

Odkładanie feromonu na krawędziach : globalne i lokalne 1) Globalne: Ma preferować najkrótszy szlak. Następuje na całej trasie gdy mrówki wykonają pełen obieg i porównają długości swoich tras. Najkrótsza w jednym obiegu: czerwona Feromon na każdej z krawędzi należącej do najkrótszej trasy zostaje zmieniony wg. f(i, j): =(1 -a)f(i, j)+a/D dla wszystkich krawędzi (i, j) na najkrótszej trasie, D – długość najkrótszej trasy a– parametr z (0, 1) a = 0 - ostatni przebieg ignorowany a=1 ignorowane doświadczenie z poprzednich przebiegów zazwyczaj optymalne około a=0. 1

Odkładanie feromonu na krawędziach : globalne i lokalne 1) Globalne (zachęcające, wykonywane po każdym obiegu) f(i, j): =(1 -a)f(i, j)+a/D D – długość najkrótszej trasy 2) Lokalne (zniechęcające, wykonywane po każdym kroku) Ma zniechęcać mrówki w przed poruszaniem się po dokładnie tej samej trasie (w tym samym obiegu) tym samym stymulując pozostałe do poszukiwań. Następuje na krawędzi (i, j) po przejściu mrówki z i do j. Feromon na krawędzi (i, j): f(i, j): =(1 -a)f(i, j)+a/p trasy p – parametr zniechęcający (p=ML, M-liczba miast, L-długość wg. heurystyki najbliższego sąsiada

Algorytm kolonii mrówek dla problemu komiwojażera wprowadzamy problem – graf pełny o N wierzchołkach. z długościami krawędzi d(i, j). Na krawędzi (i, j) ustawiamy początkowy poziom feromonu f(i, j)=1/[d(i, j)]b rozmieszczamy m mrówek w losowo wybranych miastach każda z mrówek przechodzi z miasta i do miasta j, którego jeszcze nie odwiedziła z pstwem q 0 miasto j jest wybierane tak aby t(i, j) max z pstwem 1 -q 0 miasto wybierane jest losowo z rozkładem danym funkcją preferencji. przy przejściu mrówka odkłada feromon lokalnie (zniechęcająco) na krawędzi (i, j) powtarzamy N razy każda z mrówek kończy obieg po pełnym obiegu mrówki porównują swoje trasy. feromon jest odkładany globalnie po najkrótszej trasie (zachęcająco na wszystkich krawędziach trasy) powtarzamy dopóki długość trasy ulega skróceniu

Wyniki Dobór parametru q 0 (rozkład prawdopodobieństwa zmienia się z czasem) Losowe tylko rozmieszczanie mrówek przed każdym obiegiem Eksploracja przydaje się żeby przyspieszyć zbieżność. okres nauki

Społeczność uczy się najkrótszej trasy cd: 20 miast PL południowo-wschodnia trasa optymalna – na niebiesko czerwone – krawędzie najbardziej preferowane [ max t(i, j) ] 10 iteracja Tylko 3 najbardziej preferowane krawędzie należą do trasy optymalnej. Przypadkowe zagęszczenie w okolicach Lublina i Świdnika 25 iteracja 500 iteracja Tylko 2 najbardziej preferowane krawędzie należą do trasy optymalnej: Przemyśl-Rzeszów, Katowice-Częstochowa Wszystkie najbardziej preferowane należą do trasy optymalnej

Optymalna liczba mrówek Okres nauki

Bartłomiej Urbaniec: 20 krawędzi o największej funkcji preferencji w kolejnych iteracjach 20 najczęściej używanych krawędzi parametry b=2, q 0=0. 2, m=10

Trasy mrówek w ostatniej iteracji (Bartłomiej Urbaniec)

Bartłomiej Urbaniec: Zastosowanie algorytmu dla problemu najkrótszej drogi (do problemu rozwiązywanego przez prawdziwe mrowisko – który stanowił inspirację dla stworzenia algorytmu) 1) 2) wprowadzamy prostokątna siatkę punktów punkty nie są łączone w graf pełny łączymy krawędziami tylko najbliższych sąsiadów punkt czerwony przeszkoda długość czerwonych krawędzi ustawiany nieskończona pozostałe d(i, j)=1 3) mrówka wybiera spośród (zazwyczaj) 4 możliwości ruchu 4) pozwalamy jej wracać do punktów już odwiedzonych 5) nakładamy ograniczenie na maksymalną liczbę kroków

1 obieg 50 250 750

Przewodnik do przeszukiwania losowego - inspiracje przyrodnicze Przyrodnicze (naturalne) algorytmy optymalizacji MC 1) teoria doboru naturalnego – algorytmy genetycze wygląd owada zoptymalizowany na drodze przypadkowego krzyżowania genów oraz mutacji z mechanizmem selekcji naturalnej 2) wzrost i hodowla kryształów, metalurgia – algorytm symulowanego wygrzewania

orbitale walencyjne 2 px, 2 py, 2 pz -węgiel -tworzy kierunkowe wiązania kowalencyjne atom węgla: nieobsadzone miejsca na orbitalach +6 e Stabilne formy węgla: każda – lokalne minimum energii diament E wiązania (funkcja struktury układu)=E wiązania (1023 położeń atomów) grafit E grafit diament

Wzrost kryształów jako proces optymalizacji E wiązania (funkcja struktury układu)=E wiązania (1023 położeń atomów) metoda Czochralskiego hodowli kryształów zarodek krystaliczny niska T roztopiony materiał ciut powyżej temperatury topnienia wysoka T http: //www. fkf. mpg. de

Wzrost kryształów jako proces optymalizacji E wiązania ( funkcja struktury układu)=E wiązania ( 1023 położeń atomów) metoda Czochralskiego wzrostu kryształów zarodek roztopiony materiał nieco powyżej temperatury topnienia zarodek wolno wyciągany roztopiony materiał stygnie i powoli krystalizuje jeśli odpowiednio wolno schładzany materiał krystalizuje w idealnej strukturze (o optymalnej energii wiązania) http: //www. fkf. mpg. de jeśli zarodek zbyt szybko wyciągnięty: kryształ będzie złej jakości – defekty [układ osiąga najbliższe minimum lokalne]

Struktura krystaliczna i defekty: położenie atomu wakansja dyslokacja krawędziowa energia kryształu położenie międzywęzłowe Defekty powodują naprężenia wewnętrzne. Kryształ z defektami jest twardy. przywrócenie idealnej struktury: wymaga pokonania bariery energetycznej dla usunięcia defektów (usunięcia naprężeń i zmiany twardości metalu) kryształ nagrzewa się do wysokiej temperatury, potem powoli schładza. Bariera energetyczna pokonana dzięki energii dostarczonej w formie ciepła. [proces odwrotny do hartowania stali]

Metoda z symulacji mechaniki statystycznej - układu o bardzo wielu stopniach swobody w równowadze termicznej z otoczeniem – algorytm Metropolisa. Kirkpatrick –algorytm Metropolisa do optymalizacji złożonych problemów Symulowane wygrzewanie (simulated annealing) Kirkpatrick, Science 220 671 1983 praca wykonana w IBM przy optymalizacji fizycznego projektowania układów scalonych) Pomysł: optymalizacja funkcji wielu zmiennych naśladująca proces usuwania defektów z kryształu. usuwanie defektów: optymalizacja energii wiązania w funkcji położeń wielkiej liczby atomów. Idealne optimum osiągane, gdy układ powoli schładzany (tak aby zachowana chwilowo równowaga termiczna).

Zachowanie układów o dużej liczbie stopni swobody w równowadze termicznej ze zbiornikiem ciepła kilka informacji o równowadze termodynamicznej i temperaturze – przyda się dla również dla symulacji transportu ciepła (źródło: F. Reiff Mechanika Statystyczna) kontakt termiczny

Zachowanie układów o dużej liczbie stopni swobody w równowadze termicznej ze zbiornikiem ciepła otoczenie układ Eu kontakt termiczny Eo

otoczenie układ Eu kontakt termiczny Eo całość izolowana = układ zamknięty

otoczenie układ Eu kontakt termiczny Eo całość izolowana = układ zamknięty Ec=Eu+Eo (zachowana)

otoczenie układ Eu Eo kontakt termiczny Ec – niezależna od czasu całość izolowana = układ zamknięty Ec=Eu+Eo (zachowana)

za F. Reiff Mechanika Statystyczna otoczenie układ Eu Eo kontakt termiczny całość izolowana = układ zamknięty Ec=Eu+Eo (zachowana) Ec – niezależna od czasu ale Ec=Eu(t)+Eo(t) – układy wymieniają energię wymiana energii = w stanie równowagi [równe temperatury zbiorników] fluktuacje (uśredniające się do zera) = poza równowagą [różne temperatury] ukierunkowany transfer ciepła dla wyrównania temperatur (doprowadzenia do równowagi)

dążenie do równowagi za F. Reiff Mechanika Statystyczna gorący otoczenie zimny układ Eu Eo całość izolowana = układ zamknięty kontakt termiczny Ec=Eu+Eo (zachowana) przekaz energii między układem „gorącym” i „zimnym” aż do osiągnięcia równowagi stany mikroskopowe gorącego i zimnego zbiornika gorący zimny (np. blok lodu) (np. garnek z wrzątkiem)

W(E) liczba stanów mikroskopowych układu (stan mikro – położenia i prędkości cząsteczek), dla której energia zgromadzona w układzie wynosi E stan mikroskopowy o minimalnej energii – jeden stanów mikro o większej energii – znacznie więcej założenie: każdy stan mikroskopowy jest równie prawdopodobny wtedy: pstwo, że realizowany dany stan makroskopowy = proporcjonalne do liczby stanów mikro odpowiadających danemu stanowi makro P(Eu)=CW(Eu)W’(Ec-Eu) mniejszy podukład większy W(E) – funkcja silnie rosnąca z E (oraz z rozmiarem układu) pstwo, że układ zawiera energię Eu (że podział energii Eu, Eo=Ec-Eu) (inaczej: że układ w takim stanie makroskopowym, że u ma energię Eu) równowaga (jedna z definicji): stan najbardziej prawdopodobny

P(Eu)=CW(Eu)W’(Ec-Eu) W(E) – funkcja silnie rosnąca z E (oraz z rozmiarem układu) W E u E transfer energii (w formie ciepła) między u a o aż do osiągnięcia równowagi stan równowagi: maksymalne prawdopodobieństwo. Możliwe fluktuacje wokół wartości najbardziej prawdopodobnej.

P(Eu)=CW(Eu)W’(Ec-Eu) W(E) – funkcja silnie rosnąca z E (oraz z rozmiarem układu) W jeśli P max to również log P max: poniżej na tej folii log = ln Eo=Ec- Eu w równowadze (stan najbardziej prawdopodobny): E u E transfer energii (w formie ciepła) w równowadze: między u a o aż do osiągnięcia równowagi stan równowagi: maksymalne prawdopodobieństwo. Możliwe fluktuacje wokół wartości najbardziej prawdopodobnej. (def. temperatury) transfer ciepła – napędzany gradientem T i dążący do jej wyrównania

w równowadze: układ bardzo mały w porównaniu z otoczeniem, z którym znajduje się w równowadze termicznej pytanie: jakie jest prawdopodobieństwo, że układ ma energię Eu : P(Eu) ? Eu<<Ec, W<<W’ Taylor: inne czynnik Boltzmanna

granica wysokiego T: układ może znaleźć się z równym prawdopodobieństwem przechowywać dowolnie wysoką energię 1. 0 N 0. 8 exp(-E/10) 0. 6 0. 4 0. 2 0. 0 0 exp(-E) exp(-E*10) 2 4 E 6 8 10 znalezienie układu w stanie o dużej energii jest prawdopodobne niska T- układ będzie przechowywał minimalną wartość energii jeśli T bliskie zeru: znajdziemy go tylko w stanie o najmniejszej energii

Algorytm Metropolisa – symulacja układu, który fluktuuje termicznie w równowadze z otoczeniem (zbiornikiem ciepła) kolejne kroki czasowe – zmiany stanu mikroskopowego układu układ ma podlegać rozkładowi Boltzmana. 1. 0 N 0. 8 exp(-E/10) 0. 6 0. 4 0. 2 0. 0 0 exp(-E) exp(-E*10) 2 4 E 6 8 10

Algorytm Metropolisa za S. Kooninem Układ opisany zbiorem zmiennych X=(x 1, x 2, x 3, . . . , x. N) – punkt w wielowymiarowej przestrzeni chcemy zasymulować fluktuacje termiczne zmiennych, X 1 X 2 X 3 X 4 X 5 tak, żeby układ był w „równowadze” termicznej z otoczeniem o temperaturze T. prawdopodobieństwo tego, że na „ścieżce” fluktuacji pojawi się stan X ma być dane przez zadany z góry rozkład w(X), u nas: Boltzmana w(X)= Cexp(-E(X)/k. T)

Algorytm Metropolisa X Xp układ jest w stanie X. Do jakiego stanu przejdzie? wybieramy losowo stan próbny Xp z pewnego otoczenia stanu X Xp: =X+(dx 1, dx 2, . . . , dxn) – dxi losowane z przedziału [-d, d] z rozkładem równomiernym możliwości: 1) Xp mniej prawdopodobne 2) Xp bardziej prawdopodobne w(X) [rozkład pstwa] X Xp punkt mniej prawdopodobny akceptujemy z pstwem w(Xp)/w(X) liczymy r=w(Xp)/w(X) jeśli r 1 [ w(Xp) w(X) ] – akceptujemy stan próbny X: =Xp w(X) Xp X zawsze akceptujemy bardziej prawdopodobny punkt jeśli r<1 [w(Xp)<w(X) ] akceptujemy punkt próbny z prawdopodobieństwem r losujemy liczbę losową q z przedziału [0, 1] jeśli q < r – akceptujemy stan próbny X: =Xp jeśli q r – odrzucamy

punkty na ścieżce wygenerowanej algorytmem M. istotnie podlegają rozkładowi pstwa w(X): wyobraźmy sobie, że mamy wielu wędrowców poruszających się zgodnie z algorytmem Metropolisa X Y uwagę koncentrujemy na 2 stanach układu N(X) – aktualna liczba wędrowców w punkcie X zmiana liczby wędrowców w Y w wyniku migracji z i do punktu X D = N(X)P(X Y) - N(Y)P(Y X) w warunkach równowagi (bez zmian w średnim rozkładzie wędrowców Nr(X) ) równowaga osiągana po odpowiednio dużej liczbie kroków algorytmu

równowaga w schemacie Metropolisa N=10 000 wędrowców. rozkład w(x)=C exp(-50 x 2) niebieska linia łączy punkty dane kropki: liczba wędrowców w przedziale x-dx/2, x+dx/2 dx=0. 01 start: rozkład równomierny

Algorytm Metropolisa pozostaje pokazać, że Nr(X) jest proporcjonalne do w(X) pstwo podjęcia próby (wylosowania Y jako p. próbnego) pstwo zaakceptowania próby losujemy punkty próbne z rozkładem równomiernym mamy 2 przypadki: w(X)>w(Y) : A(Y X) = 1 A(X Y) =w(Y)/w(X)<w(Y) : A(X Y) = 1 A(Y X) =w(X)/w(Y) co chcieliśmy zobaczyć

Kirkpatrick – symulowane wygrzewanie w[E(X)]=Cexp(-E(X)/k. T) alg. Metropolisa z modyfikacją rozkładu pstwa (T) w miarę działania algorytmu Algorytm symulowanego wyżarzania dla optymalizacji E(P) Wystartuj w punkcie P, ustaw wysoką „temperaturę” T Przesuń P losowo P’=P+d. P Nowy punkt akceptowany (P: =P’) zawsze gdy lepszy E(P’)<E(P) (lepszy=bardziej prawdopodobny wg. rozkładu B) gdy E(P’)>E(P) –prawdopodobieństwo zaakceptowania punktu gorszego P’ dane przez np. exp(-(E(P’)-E(P)) / k. T) (rozkład Boltzmana) zmniejszyć T Koniec jeśli T=0 Losujemy liczbę losową q wg rozkładu równomiernego, jeśli q < exp (-(E(P’)-E(P)) / k. T) P: =P’ Im niższe T, tym mniej chętnie akceptujemy przesunięcia w górę na skali energii

Przykład 1: f(x)=sin(x)+x 2/1000 Sposób zmiany temperatury: T=0. 001 i 2 gdzie i spada od 100 do 1 w każdej T wykonywana pewna liczba losowań (przesunięcia z przedziału [-2, 2]) Wysokie T – punkt P wędruje między minimami Niskie T – P uwięziony wokół jednego minimum W wysokiej T przeszukiwany szeroki zakres zmiennych. W niższej – algorytm bada dokładnie minimum, które może być globalne, jeśli schładzanie zostało odpowiednio wykonane. Techniczna uwaga: gdy zmieniamy T: najlepiej startować od najlepszego rozwiązania uzyskanego do tej pory.

Przykład 2. Zastosowanie S. A. Dla funkcji testowej de Jonga: Położenia P w kolejnych iteracjach: T>5 T<1

zawężanie zakresu przeszukiwań z temperaturą: generowane ścieżki ścieżka (100 kroków) dla T=4 ścieżka (50 kroków) dla T=2 ścieżka (25 kroków) dla T=1

Przykład 3: Problem komiwojażera Generowanie P’ z P: Losowo zmieniamy kolejność P=[13645972] Losujemy pierwsze i ostatnie miasto P’ = [ 1 3 9 4 6 5 7 2 ] Strategia schładzania a wynik końcowy Sposób zmiany T długość P zbyt szybkie schłodzenie [do najbliższego minimum lokalnego] Najlepsze rozwiązanie

Przykład 4: klaster jonowy http: //www. physchem. co. za/Bonding/Graphics/GRD 60002. gif najprostsze przybliżenie: naładowane kule bilardowe (nie mogą się przenikać) r d dodatnio naładowane – współrzędne r 1, r 2, . . . ujemnie naładowane – współrzędne d 1, d 2, . . .

Przykład 4: klaster jonowy potencjał oddziaływania: załóżmy że promienie jonowe są równe 1 oddziaływanie jonów o różnym znaku o tym samym znaku najprostsze przybliżenie: naładowane kule bilardowe (nie mogą się przenikać) r d dodatnio naładowane – współrzędne r 1, r 2, . . . ujemnie naładowane – współrzędne d 1, d 2, . . .

Weźmy 5 jonów dodatnich, pięc ujemnych, dwa wymiary r fcja 20 zmiennych ją azu k o p ałki z r t s ry cji, a u l t u a r z sym j tempe e j c ura bniżane g i f n ko la o d i k i wyn -4. 778046 -4. 853392 d -6. 856922 -2. 895126 -7. 507789

-7. 507789 minimum globalne ominięte minima lokalne -6. 413611 -7. 130378

Przykład 5: najkrótsza droga na płaszczyźnie (kontinuum, nie graf) V(x, y) meta wzniesienia: utrudnienia na trasie można je pokonać ale to kosztuje y x przeszkody start środki przeszkód

długość z definicji: całka z 1 wzdłuż trasy: funkcja kosztu: całka liczona numerycznie rozwiązanie próbne = przesunięcie każdego z zakrętów w zakresie [-. 5, . 5]* T (zależna od T) schładzanie paraboliczne poruszamy się po odcinkach prostych: liczba zakrętów = parametr rachunku

wyniki: jeden zakręt

2 zakręty 3 zakręty !

Generacja liczb losowych o zadanym (ciągłym) rozkładzie prawdopodobieństwa algorytm Metropolisa położenia „wędrowca” jest zmienną losową podległą zadanemu rozkładowi w(x) generacja liczb losowych o zadanym rozkładzie przydatne dla symulacji MC, obliczania całek, rozwiązywanie równań różniczkowych, itd.

obliczenia całek MC 1) hit or miss 1 f(x) generator o rozkładzie równomiernym 0 losujemy x losujemy y 1 0 szacowanie błędu błąd wyznaczenia średniej: odchylenie standardowe (cecha w)

obliczenia całek MC 1) hit or miss generator o rozkładzie równomiernym wartość dokładna [1 -exp(-1)] wartości plus minus e e 1/sqrt(N) N (jedna sekwencja)

1 obliczenia całek MC 2) Metoda szacowania średniej generator o rozkładzie równomiernym hit or miss f(x) 0 losujemy x 1 0 średnia hom e średnia 1/sqrt(N)

obliczenia całek MC 3) użycie niejednorodnego rozkładu pstwa taka aby g(x) mogło być g(x) 0 gęstością prawdopodobieństwa wtedy: I = <h(x)>g(x) dla zmiennej losowej x opisanej rozkładem g(x) /N jeśli g(x) tak wybrana, aby h(x) słabo zmienna w przedziale (a, b) - mamy szanse na niski błąd uwaga: jeśli h(xi)=const – błąd będzie 0

Przykład liczenia całek z użyciem a. Metropolisa 5 -5 C (gdzie C bliskie jedynce) h(x)=sqrt(5*pi) cos(x)/C <h>=<sqrt(5 p) cos(x)> N=1000 wędrowców : amplituda fluktuacji maleje z N niebieska kreska: dokładna wartość całki (1. 1355) start – wędrowcy rozmieszczeni wg rozkładu równomiernego z przedziału [-5, 5] czerwona kreska: liczba kroków a. Metropolisa