Index of /rozprawy2/10599

Pełen tekst

(1)Akademia Górniczo-Hutnicza im. Stanisława Staszica w Krakowie. Wydział Elektrotechniki, Automatyki, Informatyki i Inżynierii Biomedycznej Katedra Elektrotechniki i Elektroenergetyki. Rozprawa Doktorska Numeryczne algorytmy tomografii rezystancji siatek rezystorów. mgr inż. Piotr Zegarmistrz. Promotor: dr hab. inż. Zbigniew Galias, prof. AGH. Kraków, 2012.

(2) Składam serdeczne podziękowania prof. dr hab. inż. Zbigniewowi Galiasowi, mojemu promotorowi, za nieocenioną pomoc i wartościowe wskazówki udzielone podczas powstawania tej pracy. Dziękuję również za motywację i wiarę w moje możliwości! Szczególne podziękowania kieruję do mojej Ani, za cierpliwość, wyrozumiałość i wsparcie! Pracę dedykuję moim Rodzicom..

(3) Spis treści Skorowidz najważniejszych oznaczeń........................................................................... 4 1. Wstęp.......................................................................................................................... 6 1.1. Cel pracy ............................................................................................................ 7 1.2. Zawartość pracy................................................................................................. 8 2. Rekonstrukcja konduktancji na podstawie pomiarów brzegowych .................. 10 2.1. Siatka rezystorów............................................................................................. 10 2.2. Sformułowanie problemu ................................................................................ 12 3. Algorytmy rekonstrukcji konduktancji................................................................ 15 3.1. Metaheurystyczne i optymalizacyjne algorytmy rekonstrukcji kondunktancji15 3.1.1. Definicja funkcji celu........................................................................... 16 3.1.2. Właściwości funkcji celu ..................................................................... 19 3.1.3. Algorytmy metaheurystyczne .............................................................. 23 Metoda Monte Carlo............................................................................ 24 Symulowane wyżarzanie ..................................................................... 25 Algorytm genetyczny........................................................................... 26 3.1.4. Algorytmy optymalizacyjne ................................................................ 28 Algorytm sympleksowy Neldera-Meada ............................................. 29 Metoda quasi-Newtona ........................................................................ 31 Sekwencyjne programowanie kwadratowe ......................................... 33 Nieliniowa metoda najmniejszych kwadratów .................................... 36 3.2. Analityczny algorytm rekonstrukcji kondunktancji ........................................ 38 3.2.1. Rekonstrukcja konduktancji w jednej warstwie .................................. 39 3.2.2. Przypadek k = 1.................................................................................... 41 3.2.3. Przypadek k większego od 1 i nie większego od liczby wierszy......... 42 Wyznaczenie współczynników α1, α2, ..., αk ........................................ 42 Wyznaczenie potencjałów ................................................................... 43 Rekonstrukcja kondunktancji .............................................................. 44 3.2.4. Przypadek k większego od liczby wierszy........................................... 44 3.3. Modyfikacje algorytmu analitycznego ............................................................ 46 3.3.1. Uruchomienie algorytmu z wszystkich narożników siatki .................. 46 3.3.2. Optymalizacja wartości współczynników α1, α2, ..., αk ....................... 48 4. Wyniki działania algorytmów rekonstrukcji ....................................................... 50 4.1. Szczegóły implementacji ................................................................................. 50.

(4) 4.2. Algorytmy metaheurystyczne .......................................................................... 51 4.3. Algorytmy optymalizacyjne ............................................................................ 58 4.4. Algorytm analityczny ...................................................................................... 66 4.5. Algorytm analityczny z usprawnieniami ......................................................... 69 4.6. Porównanie metod rekonstrukcji ..................................................................... 76 5. Podsumowanie......................................................................................................... 78 Dodatek A. Obliczenia analityczne dla przykładowej siatki rezystancyjnej........... 83 Literatura ...................................................................................................................... 92 Spis rysunków i tabel.................................................................................................... 97.

(5) Skorowidz najważniejszych oznaczeń N1. –. liczba kolumn węzłów siatki rezystancyjnej (rozmiar poziomy siatki). N2. –. liczba wierszy węzłów siatki rezystancyjnej (rozmiar pionowy siatki). Nb. –. liczba węzłów brzegowych siatki rezystancyjnej. Ne. –. liczba elementów (krawędzi) siatki rezystancyjnej. Γ Ω0. –. siatka rezystancyjna. –. zbiór węzłów siatki rezystancyjnej. Ω1. –. zbiór krawędzi siatki rezystancyjnej. Pij. –. węzeł siatki prostokątnej o współrzędnych (i,j). int Ω0. –. zbiór węzłów wewnętrznych siatki rezystancyjnej. ∂Ω0. –. zbiór węzłów brzegowych siatki rezystancyjnej. γ. –. konduktywność (odwzorowanie). Ib. –. wektor prądów w węzłach brzegowych (zestaw danych Neumanna). Vb. –. wektor potencjałów węzłów brzegowych (zestaw danych Dirichleta). Λ. –. macierz odwzorowania Dirichlet-Neumann. Λ. –. macierz odwzorowania Neumann-Dirichlet. V. –. wektor potencjałów węzłów siatki. G. –. wektor konduktancji elementów siatki, G = [G1, G2, ..., GNe]. f(G). –. pA. –. ΔE. –. funkcja celu wykorzystywana w algorytmach rekonstrukcji bazujących na procedurach metaheurystycznych i optymalizacyjnych prawdopodobieństwo akceptacji kolejnej próbki (algorytm symulowanego wyżarzania) różnica energii analizowanej i poprzedniej próbki (algorytm symulowanego wyżarzania). T. –. temperatura (algorytm symulowanego wyżarzania). T0. –. początkowa wartość temperatury (algorytm symulowanego wyżarzania). tf. –. maksymalna liczba iteracji (algorytm symulowanego wyżarzania). fp(G). –. funkcja przystosowania osobnika (algorytm genetyczny). fmax. –. wartość funkcji celu f(G) dla najgorszego osobnika w pokoleniu (algorytm genetyczny). pc. –. prawdopodobieństwo krzyżowania (algorytm genetyczny). pm. –. prawdopodobieństwo mutacji (algorytm genetyczny). h. –. maksymalny względny błąd rekonstrukcji. −1.

(6) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. 1.. Wstęp. Problematyka tomografii impedancyjnej jest przedmiotem intensywnych badań od wielu lat [13, 22, 50, 51]. Dynamiczny rozwój technik tomografii impedancyjnej, zarówno dwu- jak i trójwymiarowej jest związany z możliwością licznych aplikacji praktycznych, przede wszystkim w medycynie [1], ale również w wielu innych dziedzinach, takich jak np. budownictwo czy konserwacja zabytków [21, 35]. Istotą tomografii, w każdej z jej odmian, jest określenie parametrów fizycznych danego ośrodka wyłącznie na podstawie pomiarów na jego brzegu. Większość stosowanych obecnie technik tomograficznych dotyczy badania ośrodków ciągłych. Na podstawie pomiarów brzegowych konstruowany jest obraz właściwości elektrycznych wnętrza badanego ośrodka. Dalej, na podstawie powstałego obrazu wnioskuje się na temat właściwości fizycznych ośrodka. Poszczególne techniki tomografii impedancyjnej różnią się od siebie m.in. wyborem sygnałów elektrycznych wymuszających oraz sygnałów będących odpowiedzią badanego ośrodka na wymuszenie, metodami opisu matematycznego zjawisk fizycznych występujących w procesie tomografii oraz metodami numerycznymi rozwiązywania sformułowanych problemów. W poniższej pracy analizowany jest problem tomografii rezystancji siatek rezystorów [9, 10, 17, 28, 29]. Jest to problem dyskretny, w którym na podstawie pomiarów elektrycznych na brzegu układu składającego się ze skończonej liczby elementów określane są wartości konduktancji (bądź rezystancji) tych elementów. Problem ten ma potencjalne zastosowania w dziedzinie tomografii ośrodków ciągłych. W zastosowaniach tych wyznaczane są właściwości fizyczne ośrodka ciągłego przez zamodelowanie go odpowiednio gęstym układem dyskretnym. Realizacja taka jest możliwa pod warunkiem opracowania efektywnych technik tomografii rezystancyjnej układów dyskretnych składających się z wystarczająco dużej liczby elementów. Rozważany problem może mieć również praktyczne zastosowania w sytuacjach, gdy siatka elementów rezystancyjnych jest wykorzystywana do wyznaczenia rozkładu ciężaru, temperatury lub innych parametrów fizycznych danej powierzchni. W tych zastosowaniach tworzy się siatkę rezystorów o znanym rozmiarze i strukturze połączeń. Elementy łączące węzły siatki wykonane są z tzw. „inteligentnych” materiałów, czyli. 6.

(7) Numeryczne algorytmy tomografii rezystancji siatek rezystorów takich, których właściwości elektryczne zależą od innych parametrów fizycznych. Parametrami tymi mogą być m.in. temperatura elementu lub nacisk na jego powierzchnię. W sytuacji, gdy znana jest zależność konduktancji od zmiany interesującej nas wielkości fizycznej (np. temperatury) można zbudować układ pozwalający na monitoring zmian rozkładu tego parametru fizycznego na powierzchni układu. Możliwe jest wówczas wykorzystanie metod tomografii siatek rezystorów do stworzenia obrazu rozkładu danego parametru fizycznego i obserwowanie jego zmian w czasie rzeczywistym.. 1.1. Cel pracy W swojej pracy dydaktycznej autor niniejszej dysertacji prowadził między innymi zajęcia z przedmiotów „Teoria Obwodów Elektrycznych” oraz „Computeraided Analysis of Electronic Systems”. Podczas zajęć tych, wprowadzając metodę potencjałów węzłowych, autor postawił studentom zadanie wyznaczenia wypadkowej rezystancji pomiędzy dwoma węzłami siatki rezystorów [57, 60]. W związku z tym doświadczeniem autor zainteresował się problemem rekonstrukcji konduktancji elementów siatek rezystancyjnych i podjął próbę zweryfikowania skuteczności działania istniejących oraz opracowania nowych algorytmów rekonstrukcji. Autor postanowił zbadać, czy możliwe jest wykorzystanie tych algorytmów do rekonstrukcji konduktancji siatek rezystorów w sytuacjach praktycznych oraz w tomografii ośrodków ciągłych z satysfakcjonującą dokładnością. Na podstawie badań wykonanych podczas przygotowania niniejszej dysertacji autor sformułował poniższą tezę: Możliwe jest takie sformułowanie postawionego w pracy problemu, które pozwala na realizację algorytmów rekonstrukcji bazujących na metodach optymalizacyjnych poszukiwania ekstremum funkcji wielu zmiennych. Istniejące algorytmy rekonstrukcji konduktancji elementów skończonych siatek rezystorów na podstawie pomiarów brzegowych cechuje ograniczenie dotyczące rozmiarów siatek możliwych do rekonstrukcji. Możliwe jest ulepszenie istniejących algorytmów rekonstrukcji w taki sposób, aby pozwalały na rekonstrukcję układów składających się z większej liczby elementów oraz były bardziej odporne na występowanie niedokładności pomiarów.. 7.

(8) Numeryczne algorytmy tomografii rezystancji siatek rezystorów W celu wykazania prawdziwości postawionych tez autor wyznaczył sobie następujące cele cząstkowe: – –. –. – –. Dokonać przeglądu literaturowego istniejących metod rekonstrukcji konduktancji siatek rezystorów, zaimplementować i przetestować wybrane metody. Zdefiniować funkcję celu postawionego problemu w taki sposób, aby określała ona stopień dopasowania badanej próbki/siatki testowej do rekonstruowanego oryginału. Funkcja ta powinna być tak skonstruowana, aby jej argumentem był wektor konduktancji wszystkich elementów siatki i funkcja przyjmowała ekstremum globalne dla przypadku, gdy wektor konduktancji badanej próbki jest tożsamy z wektorem konduktancji oryginalnej siatki. Na podstawie zdefiniowanej funkcji celu zaimplementować algorytmy rekonstrukcji bazujące na metaheurystycznych oraz optymalizacyjnych metodach poszukiwania minimum funkcji. Zdefiniować parametry porównawcze pozwalające ocenić skuteczność i precyzję działania poszczególnych algorytmów rekonstrukcji konduktancji. Zaimplementować i przetestować zaproponowane algorytmy, porównać ich właściwości z istniejącymi metodami rekonstrukcji.. –. Wskazać możliwe modyfikacje zaimplementowanych algorytmów pozwalające na poprawienie ich skuteczności. Oczekiwane rezultaty zweryfikować badaniami symulacyjnymi. – Zweryfikować możliwość wykorzystania algorytmów rekonstrukcji konduktancji siatek rezystorów w zastosowaniach praktycznych oraz do tomografii ośrodków ciągłych.. 1.2. Zawartość pracy Praca składa się z pięciu rozdziałów, dodatku A, skorowidza najważniejszych oznaczeń oraz spisu literatury. Pierwszy rozdział stanowi wstęp, opisano w nim pojęcie tomografii i przedstawiono możliwe zastosowania praktyczne. Ponadto przedstawiono genezę i cel pracy, sformułowano tezę pracy, zdefiniowano postawione cele cząstkowe. Streszczono także zawartość pracy. W rozdziale drugim sformułowano i szczegółowo opisano postawiony w pracy problem. Zdefiniowano podstawowe parametry skończonych, prostokątnych siatek rezystancyjnych. Pokrótce przedstawiono tematykę problemów odwrotnych. Ponadto. 8.

(9) Numeryczne algorytmy tomografii rezystancji siatek rezystorów wskazano odwzorowania w jednoznaczny sposób opisujące siatkę rezystancyjną oraz właściwości takich odwzorowań. Rozdział trzeci stanowi opis algorytmów rekonstrukcji konduktancji dla siatek prostokątnych. Zdefiniowano funkcję celu umożliwiającą przedstawienie problemu rekonstrukcji jako problemu poszukiwania minimum funkcji wielu zmiennych. Wskazano właściwości omawianej funkcji celu. Dalej zaprezentowano wybrane metody metaheurystyczne i optymalizacyjne wykorzystywane w zaproponowanych algorytmach. Następnie przedstawiono podstawowe założenia oryginalnego algorytmu analitycznego Curtisa i Morrowa oraz jego rozszerzenia pozwalające na rekonstrukcję siatek prostokątnych. W końcowej części rozdziału trzeciego przedstawiono modyfikacje algorytmu analitycznego pozwalające na zwiększenie zakresu rozmiarów siatek możliwych do rekonstrukcji oraz zmniejszenie wpływu błędów pomiarowych na wynik działania algorytmu. W rozdziale czwartym zebrano wyniki przeprowadzonych badań symulacyjnych. Zaprezentowano analizę skuteczności działania poszczególnych metod. Parametrami podlegającymi ocenie są uśredniony maksymalny błąd rekonstrukcji oraz całkowity czas wykonywania algorytmu. W przypadku algorytmów analitycznych bazujących na idei oryginalnego algorytmu Curtisa i Morrowa przedstawiono wyniki badania odporności algorytmów na występowanie błędów pomiarowych oraz wpływ niedokładności pomiaru na maksymalny rozmiar siatek możliwych do rekonstrukcji. Podsumowanie niniejszej dysertacji zawarto w rozdziale piątym. Zostały tam zaprezentowane wnioski wyciągnięte na podstawie przeprowadzonych badań. Wskazane zostały możliwe aplikacje praktyczne opracowanych i przebadanych algorytmów. Wymienione zostały zadania wykonane w trakcie powstawania niniejszej pracy. W ostatniej części autor podsumował realizację celów cząstkowych i ustosunkował się do postawionej we wstępie tezy pracy. W Dodatku A zaprezentowano obliczenia analityczne dla przykładowej siatki o rozmiarze 4×4 węzły, składającej się z 12 elementów. Obliczono, jak współczynniki macierzy odwzorowania Neumann – Dirichlet zależą od konduktancji elementów siatki. Ponadto zobrazowano, jak wartość funkcji celu zdefiniowanej w rozdziale trzecim oraz jej pochodnych cząstkowych zależy od konduktancji elementów siatki.. 9.

(10) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. 2.. Rekonstrukcja konduktancji na podstawie pomiarów brzegowych. W pracy analizowane są algorytmy rekonstrukcji konduktancji elementów łączących węzły skończonych, prostokątnych siatek rezystancyjnych. W niniejszym rozdziale szczegółowo zdefiniowano siatki rezystancyjne oraz sformułowano postawiony w pracy problem.. 2.1. Siatka rezystorów Siatkę rezystorów można zdefiniować jako obwód elektryczny, w którym elementy rezystancyjne łączą sąsiednie węzły kratownicy. W pracy rozważany jest przypadek, gdy węzły siatki wypełniają kratownicę o kształcie prostokąta. Założono, że każdy węzeł wewnętrzny jest połączony z czterema węzłami sąsiednimi za pomocą elementu o zadanej konduktancji, która jest skończona i dodatnia, podczas gdy każdy węzeł brzegowy jest połączony z tylko jednym węzłem wewnętrznym. Obrazuje to przykład przedstawiony na rys. 2.1.. Rys. 2.1. Siatka rezystorów o wymiarach N1×N2=8×5 węzłów.. 10.

(11) Numeryczne algorytmy tomografii rezystancji siatek rezystorów Przyjmując, że węzły siatki są rozmieszczone w N1 kolumnach i N2 wierszach, można zauważyć, że siatka taka posiada (N1–2)(N2–2) węzłów wewnętrznych oraz 2(N1+N2–4) węzłów brzegowych, co daje całkowitą liczbę N1N2–4 węzłów. Dla przedstawionego przykładu siatki o rozmiarze 8×5 węzłów oznacza to 18 węzłów wewnętrznych oraz 18 węzłów brzegowych. W pracach [9, 10] przedstawiono matematyczną definicję kwadratowych siatek rezystancyjnych. Na tej podstawie poniżej przedstawiono pełną definicję prostokątnej siatki rezystancyjnej. Prostokątną siatkę rezystancyjną Γ można zdefiniować jako zbiór węzłów Ω0 , krawędzi łączących węzły Ω1 wraz z odwzorowaniem γ , tzn.: w2.1. Γ = (Ω0 , Ω1 , γ ). (2.1). γ : Ω1 → R +. (2.2). przy czym: w2.2. jest odwzorowaniem przyporządkowującym każdej krawędzi liczbę rzeczywistą dodatnią. Węzłami siatki nazywamy zbiór punktów Pij prostokątnej kratownicy o współrzędnych całkowitych, takich że Pij = (i,j), a 1 ≤ i ≤ N1 oraz 1 ≤ j ≤ N2, z wyłączeniem punktów narożnych (1,1), (N1,1), (1,N2) oraz (N1,N2). Zbiór węzłów wewnętrznych siatki, oznaczony int Ω0 , zawiera te węzły, dla których 2 ≤ i ≤ N1–1 oraz 2 ≤ j ≤ N2–1. Zbiór węzłów brzegowych siatki, oznaczony ∂Ω0 , to zbiór zdefiniowany jako Ω0 – int Ω0 . Każdy węzeł wewnętrzny ma dokładnie cztery węzły sąsiednie, przy czym odległość danego węzła od dowolnego z węzłów sąsiednich wynosi dokładnie jeden. Każdy węzeł brzegowy ma tylko jeden węzeł sąsiedni i jest to węzeł oddalony od niego dokładnie o jeden należący do wnętrza siatki. Krawędzią PijPkl nazywane jest połączenie pary węzłów Pij oraz Pkl, gdzie |k–i|+|l–j|=1. Zbiór krawędzi siatki rezystancyjnej oznaczony jest symbolem Ω1 . Z własności węzła brzegowego wynika, że nie istnieją krawędzie łączące węzły brzegowe. Dla każdej krawędzi σ ∈ Ω1 wartość funkcji γ (σ ) jest nazywana konduktancją elementu, a jej odwrotność – rezystancją. Odwzorowanie γ przyporządkowujące konduktancję każdej krawędzi kratownicy nazywane jest konduktywnością. Należy zwrócić uwagę, że kształt siatki rezystorów przedstawiony na rys. 2.1 nie jest ograniczeniem. W rozważanym przypadku węzły brzegowe nie są połączone ze sobą. Jeżeli przyjąć założenie odwrotne, tzn. siatka zawiera połączenia o niezerowej konduktancji pomiędzy węzłami brzegowymi, można wtedy dołączyć do wszystkich węzłów brzegowych dodatkowe konduktancje o ustalonej wartości. Wtedy siatka uzyska omawiany kształt, a liczby wierszy i kolumn zwiększą się o 2.. 11.

(12) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. 2.2. Sformułowanie problemu W analizie teoretycznej obwodów elektrycznych mamy z reguły do czynienia z problemami postawionymi wprost, tzn. przedmiotem analizy jest wyznaczenie odpowiedzi układu o znanej topologii i znanych wartościach elementów na zadane wymuszenie. Zagadnienie opisane w niniejszej pracy zalicza się do grupy tzw. problemów odwrotnych (ang. inverse problem). Istotą tego typu zagadnień jest wyznaczenie parametrów elementów układu (w przypadku układów dyskretnych) lub własności fizycznych obiektu (w przypadku układów ciągłych) na podstawie odpowiedzi układu na pewne wymuszenie lub zestaw wymuszeń. Postawiony w niniejszej pracy problem odwrotny polega na wyznaczeniu wartości konduktancji wszystkich elementów prostokątnej siatki rezystancyjnej na podstawie pomiarów brzegowych. Zakłada się, że istnieje fizyczny dostęp jedynie do węzłów należących do brzegu siatki ∂Ω0 . Można sobie wyobrazić dwie metody uzyskania informacji na temat konduktancji elementów siatki rezystancyjnej. W pierwszej wersji, wymuszeniem jest wektor potencjałów węzłów brzegowych. Odpowiedzią układu na to wymuszenie jest wektor prądów w elementach brzegowych. Elementem brzegowym nazywamy element łączący węzeł brzegowy z sąsiadującym z nim węzłem wewnętrznym. W drugiej wersji mamy do czynienia z sytuacją odwrotną. Wektor potencjałów w węzłach brzegowych jest odpowiedzią układu na wymuszenie w postaci wektora prądów w elementach brzegowych. W literaturze wektor prądów brzegowych siatki rezystancyjnej zwany jest zestawem danych Neumanna, podczas gdy wektor potencjałów węzłów brzegowych zwany jest zestawem danych Dirichleta. Dla każdej prostokątnej siatki rezystancyjnej można wyznaczyć współczynniki macierzy Λ będącej odwzorowaniem Dirichlet-Neumann, czyli wyznaczającej odpowiedź prądową układu Ib na wymuszenie w postaci wektora potencjałów brzegowych Vb, tj.: w2.3. Ib = Λ ·Vb. (2.3). Λ−1 będącą odwzorowaniem Neumann-Dirichlet, dzięki której można wyznaczyć odpowiedź Vb układu na wymuszenie wektorem prądów brzegowych Ib. Analogicznie, dla każdej siatki można zapisać macierz. w2.4.. Vb = Λ−1 ·Ib. (2.4). Siatka rezystancyjna o N1 kolumnach, N2 wierszach węzłów i kształcie przestawionym na rys. 2.1 ma Nb = 2(N1+N2–4) węzłów brzegowych. Macierze Λ oraz Λ−1 są rozmiaru (Nb–1)×(Nb–1), czyli składają się z (Nb–1)2 elementów. Wymiar macierzy Λ oraz Λ−1 jest o jeden mniejszy od liczby węzłów brzegowych. Wynika to. 12.

(13) Numeryczne algorytmy tomografii rezystancji siatek rezystorów bezpośrednio z wymiaru wektorów Vb oraz Ib. Wymiar wektora Vb wynosi Nb –1, gdyż zawsze przyjmuje się jeden z węzłów brzegowych, jako węzeł odniesienia o zerowym potencjale. Potencjały pozostałych węzłów określane są względem węzła odniesienia. Wymiar wektora Ib również wynosi Nb –1. Jeżeli wymusimy prądy we węzłach brzegowych poza jednym, to wartość prądu w tym ostatnim węźle wynika z prądowego prawa Kirchhoffa. Dla przykładu z rys. 2.1 macierze Λ oraz Λ−1 składają się z 172 = 289 elementów. W pracach [9, 10] udowodniono, że problem rekonstrukcji ma jednoznaczne rozwiązanie, tzn. że odwzorowania Dirichlet-Neumann oraz Neumann-Dirichlet w sposób jednoznaczny określają wartości wszystkich elementów siatki rezystancyjnej. Ponadto wykazano, że problem ten ma charakter ciągły, tzn. że jeżeli wartości pomiarów na brzegach siatki są do siebie zbliżone, to siatki, dla których uzyskano te wyniki są zbudowane z elementów o zbliżonych wartościach konduktancji. Elementy macierzy Λ lub Λ−1 można wyznaczyć na podstawie skończonej liczby pomiarów. Przykładowo, możliwe jest wyznaczenie wszystkich współczynników macierzy Λ na podstawie Nb–1 zestawów pomiarowych. W każdym kroku rozpływ prądów w siatce jest wymuszony wektorem potencjałów brzegowych Vb, w którym tylko jeden element jest różny od zera, tzn. przykładamy napięcie E = 1V do jednego z węzłów brzegowych, pozostałe Nb–1 węzłów uziemiamy i dokonujemy pomiaru prądów brzegowych, zapisując je w wektorze Ib. W ten sposób w i-tym kroku tak opisanego algorytmu wyznaczane są wartości współczynników i-tej kolumny macierzy Λ . Po Nb–1 krokach znane więc są wartości wszystkich elementów macierzy odwzorowania Dirichlet-Neumann Λ . Współczynniki macierzy odwzorowania Neumann-Dirichlet Λ−1 można wyznaczyć w podobny sposób. W kolejnych krokach należy wymuszać rozpływ prądów w obwodzie wektorami Ib, w których występuje tylko jeden niezerowy element. W wybranym węźle brzegowym jest wymuszany prąd o wartości J = 1A, podczas gdy w pozostałych węzłach brzegowych, poza węzłem odniesienia, prąd wynosi zero. Procedura wymaga wykonania Nb–1 pomiarów, przy czym przy każdym pomiarze niezerowy prąd jest wymuszany w innym węźle brzegowym. Wyznaczenie węzła odniesienia gwarantuje spełnienie prądowego prawa Kirchhoffa przy wymuszaniu rozpływu prądów w siatce opisanymi powyżej wektorami Ib. Ponadto warto zauważyć, że zgodnie z wzorami (2.3) i (2.4) macierz odwzorowania Neumann-Dirichlet Λ−1 jest macierzą odwrotną macierzy odwzorowania Dirichlet-Neumann Λ . Postawiony w pracy problem można przedstawić jako zagadnienie wyznaczenia wartości konduktancji wszystkich elementów siatki rezystancyjnej na podstawie. 13.

(14) Numeryczne algorytmy tomografii rezystancji siatek rezystorów współczynników macierzy Λ lub Λ−1 . Liczba elementów siatki w funkcji jej rozmiaru wynosi: w2.5.. Ne=(N1–1)(N2–2)+(N1–2)(N2–1)=2N1N2–3N1–3N2+4. (2.5). natomiast liczba elementów macierzy Λ lub Λ−1 wynosi: w2.6.. (Nb–1)2 = (2(N1+N2–4)–1)2. (2.6). Zatem, dla każdej siatki liczba elementów macierzy Λ lub Λ−1 jest większa od liczby elementów siatki. Dowód tego, że informacja zawarta we współczynnikach macierzy Λ (lub Λ−1 ) jest wystarczająca do wyznaczenia wartości konduktancji wszystkich elementów siatki rezystancyjnej znajduje się w pracy [9]. W Dodatku A do niniejszej rozprawy przedstawiono obliczenia analityczne dla przykładowej siatki o rozmiarze 4×4 węzły. Wyprowadzono tam m.in. zależność wszystkich współczynników macierzy Λ−1 od konduktancji elementów łączących węzły siatki. Zależność ta jest bardzo skomplikowana już dla niewielkich rozmiarów siatek.. 14.

(15) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. 3.. Algorytmy rekonstrukcji konduktancji. W niniejszym rozdziale zostały szczegółowo omówione wszystkie zaimplementowane algorytmy rekonstrukcji konduktancji siatek rezystorów. W pierwszej części rozdziału zaprezentowano grupę algorytmów metaheurystycznych i optymalizacyjnych, których istota polega na poszukiwaniu najlepszego rozwiązania postawionego w pracy problemu na podstawie odpowiednio zdefiniowanej funkcji celu. Następnie przedstawiono algorytm analityczny, pozwalający na dokładne wyliczenie wartości konduktancji wszystkich elementów siatki. Algorytm ten składa się ze skończonej liczby kroków, w których konduktancje są rozpoznawane odpowiednimi warstwami. W ostatniej części rozdziału opisano zaproponowane przez autora modyfikacje algorytmu analitycznego.. 3.1. Metaheurystyczne i optymalizacyjne algorytmy rekonstrukcji kondunktancji W tej części pracy opisany został sposób wykorzystania metod metaheurystycznych i optymalizacyjnych w algorytmach rekonstrukcji konduktancji siatek rezystorów. Wspomniane metody zostały opisane razem, gdyż mają one bardzo istotną wspólną własność. W celu zrealizowania w oparciu o nie algorytmów rekonstrukcji konduktancji należy zdefiniować funkcję celu, która pozwoli na ocenę dopasowania wybranej próbki (siatki testowej) do oryginału. Funkcja ta musi być tak zdefiniowana, aby jej wartość zależała od wartości konduktancji wszystkich elementów siatki poddawanej ocenie i była miarą jej dopasowania do oryginału. Istotą algorytmów metaheurystycznych jest przeszukiwanie zbioru możliwych rozwiązań i obliczanie wartości funkcji celu dla każdej badanej próbki. W przypadku algorytmów działających w oparciu o analizę pojedynczej próbki, w każdym kroku, na podstawie wartości funkcji celu, podejmowana jest decyzja, czy analizowana próbka jest lepiej dopasowana do oryginalnej siatki od dotychczasowego najlepszego rozwiązania. Jeżeli tak jest, badana siatka testowa jest zapisywana jako nowe najlepsze. 15.

(16) Numeryczne algorytmy tomografii rezystancji siatek rezystorów rozwiązanie. W przypadku algorytmów analizujących wiele próbek w każdym kroku, na podstawie wartości funkcji celu wszystkich próbek, wybierana jest reprezentacja, która stanowić będzie bazę grupy analizowanej w kolejnym kroku. Poszczególne algorytmy metaheurystyczne różnią się metodą zmierzania do globalnego minimum. Często algorytmy metaheurystyczne inspirowane są zjawiskami fizycznymi lub przyrodniczymi. W przypadku algorytmów bazujących na metodach optymalizacyjnych należy sformułować problem rekonstrukcji konduktancji jako problem poszukiwania minimum globalnego funkcji wielu zmiennych. W tym celu należy zdefiniować funkcję celu, której zmiennymi są wartości konduktancji wszystkich elementów siatki, a jej minimum globalne odpowiada rozwiązaniu problemu rekonstrukcji. Następnie, stosując wybraną metodę optymalizacyjną, znaleźć wektor konduktancji elementów siatki dla którego funkcja celu przyjmuje wartość minimalną. Z powyższego wynika, że funkcję celu można zdefiniować w taki sposób, aby spełniała założenia realizacji algorytmów rekonstrukcji konduktancji elementów siatek rezystancyjnych, działających zarówno w oparciu o metody metaheurystyczne jak i optymalizacyjne.. 3.1.1. Definicja funkcji celu. Funkcję celu zdefiniowano w oparciu o odwzorowanie Neumann-Dirichlet, czyli odpowiedź układu w postaci wektora potencjałów węzłów brzegowych na wymuszenie w postaci wektora prądów w elementach brzegowych. W każdej serii pomiarów zakładamy, że wymuszamy odpowiednie prądy w elementach brzegowych (dane wejściowe – zestaw Neumanna) i dokonujemy pomiaru napięć pomiędzy węzłami brzegowymi (dane wyjściowe – zestaw Dirichleta). Liczba danych w wektorze potencjałów brzegowych (równa Nb–1, gdzie Nb jest liczbą węzłów brzegowych) jest zawsze mniejsza niż liczba elementów siatki rezystancyjnej. Oczywistym zatem jest, że należy wykonać więcej niż jeden zestaw pomiarów, aby zapewnić jednoznaczność rozwiązania. Ponadto należy wybrać takie wektory prądów brzegowych, aby zagwarantować, że prąd w każdym elemencie brzegowym jest różny od zera dla co najmniej jednego wektora. Jeżeli w żadnym z wybranych wektorów prądów brzegowych nie będzie wymuszony niezerowy prąd w wybranym elemencie brzegowym Gm (innymi słowy ani raz nie popłynie przez ten element prąd) to prawidłowa rekonstrukcja wartości konduktancji tego elementu nie będzie możliwa.. 16.

(17) Numeryczne algorytmy tomografii rezystancji siatek rezystorów Funkcję celu zdefiniowano w następujący sposób: l. wy1. f (G ) =. Nb − 1. ∑ ∑ (V s=1. org s,t. t=1 Nb − 1 l. ∑ ∑ (V s=1. − V tests,t )2 (3.1) org 2 s,t. ). t=1. gdzie G = [G1, G2, ..., GNe] jest wektorem konduktancji elementów siatki, Ne jest liczbą elementów siatki, l jest liczbą zestawów pomiarowych, Nb–1 jest rozmiarem wektora potencjałów brzegowych, V orgs,t oznacza zmierzoną wartość potencjału węzła t w s-tym zestawie pomiarów oryginalnej, nieznanej siatki, a V tests,t jest obliczoną wartością potencjału węzła t w s-tym zestawie pomiarów testowej siatki rezystancyjnej opisanej wektorem konduktancji G.. Rys. 3.1. Zestaw Neumanna dla pierwszego kroku obliczenia wartości funkcji celu.. W ostatecznej wersji algorytmu, do obliczenia wartości funkcji celu dla wybranej siatki testowej użyto l=N1+N2−4 zestawów pomiarowych, gdzie N1 i N2 są odpowiednio liczbami kolumn i wierszy siatki. W każdym zestawie pomiarowym tylko dwa elementy wektora wymuszeń prądowych mają niezerowe wartości. Zadany prąd (przykładowo J=1) jest wstrzykiwany do wybranego węzła brzegowego. 17.

(18) Numeryczne algorytmy tomografii rezystancji siatek rezystorów i pobierany z węzła leżącego symetrycznie po przeciwnej stronie siatki. Zestaw prądów brzegowych dla pierwszego zestawu pomiarowego zaprezentowano na rys. 3.1. W następnym kroku źródło prądowe z górnego wiersza przesuwane jest o jedną pozycję w prawo, a źródło w dolnym wierszu o jedną pozycję w lewo – symetrycznie względem środka siatki. Po osiągnięciu ostatniego węzła w górnym i dolnym wierszu (krok N1−2), w kolejnym kroku należy przesunąć źródła do górnego węzła skrajnej prawej kolumny i dolnego węzła skrajnej lewej kolumny. Ustawienie źródeł w N1−1-szym kroku procedury obliczania wartości funkcji celu przedstawiono na rys. 3.2. Ostatecznie, po osiągnięciu dolnego węzła prawej kolumny i górnego wiersza lewej kolumny (ostatni krok) obliczenie wartości funkcji celu jest zakończone.. Rys. 3.2. Zestaw Neumanna dla N1−1-szego kroku obliczenia wartości funkcji celu.. Testowano kilka wersji algorytmów wyboru wektorów wymuszeń prądowych w węzłach brzegowych przy zachowaniu wzoru (3.1) do obliczenia wartości funkcji celu dla wybranej siatki testowej. Do ostatecznej analizy wybrano jednak wersję opisaną powyżej. Przykładowe wyniki dla innych przebadanych wersji wraz z opisem algorytmu wyboru wektorów wymuszeń prądowych opisano w rozdziale 4. Gdy siatka testowa jest taka sama jak siatka oryginalna, funkcja celu (3.1) przyjmuje zerową wartość. Gdy siatka testowa różni się od siatki oryginalnej, różnią się również wektory potencjałów brzegowych porównywanych siatek, będące odpowiedziami na ten sam wektor prądów brzegowych. Funkcja celu przyjmuje więc wartość dodatnią. 18.

(19) Numeryczne algorytmy tomografii rezystancji siatek rezystorów Taka definicja funkcji celu pozwala na miarodajną ocenę każdej analizowanej siatki testowej. Ponadto, problem rekonstrukcji konduktancji siatki można sformułować jako problem znalezienia minimum funkcji wielu zmiennych f.. 3.1.2. Właściwości funkcji celu. Jak wspomniano powyżej, funkcja celu jest funkcją wielu zmiennych. Liczba zmiennych jest równa liczbie elementów siatki. Oznacza to, że już dla niewielkich rozmiarów siatek mamy do czynienia z problemem wysokowymiarowym. W związku z tym niemożliwa jest bezpośrednia wizualizacja zależności wartości funkcji celu od wszystkich parametrów wejściowych. Aby określić pewne charakterystyczne właściwości funkcji celu wygenerowano wykresy przedstawiające zależność wartości funkcji celu od dwóch wybranych konduktancji siatki. Wartość funkcji celu jest obliczana dla siatek, w których wszystkie konduktancje poza dwiema wybranymi do analizy mają wartości równe wartościom konduktancji oryginalnej siatki, a wybrane dwa parametry są zmieniane w całym zakresie wybranego przedziału z zadaną rozdzielczością. Punkt środkowy płaszczyzny X-Y odpowiada sytuacji, gdy wartości wybranych konduktancji są równe wartościom dla oryginalnej siatki. W punkcie tym funkcja celu przyjmuje wartość zero. W pozostałych punktach funkcja celu przyjmuje wartości dodatnie. Przykładowe wykresy przedstawiono na rys. 3.3 – 3.10. Najważniejszą obserwacją jest, że funkcja celu może być bardzo płaska w otoczeniu minimum. Wynika z tego, że znalezienie w przestrzeni poszukiwań punktu o bardzo małej wartości funkcji celu nie oznacza, że znalezione rozwiązanie jest blisko oryginału. Inną istotną obserwacją jest, że kształt uzyskanych wykresów bardzo mocno zależy od wyboru zmiennych, od których uzależniamy wartość funkcji celu. Niektóre konduktancje siatki mają duży wpływ na wartość funkcji celu, podczas gdy w przypadku innych wpływ ten jest praktycznie niezauważalny. Można zauważyć, że wartości konduktancji leżących w pobliżu brzegu siatki mają duży wpływ na wartość funkcji celu, podczas gdy wartości konduktancji leżących w pobliżu środka nie wpływają na nią znacząco. Oznacza to, iż poprawność rekonstrukcji konduktancji zależy od położenia elementów w siatce. Im bliżej środka siatki położony jest dany element, tym trudniej jest wyznaczyć jego konduktancję z satysfakcjonującą dokładnością.. 19.

(20) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. 0.8 0.7 0.6. f(G). 0.5 0.4 0.3 0.2 0.1 0 10 8. 10 6. 8. G17. 6. 4. 4. 2. G. 10. 2 0. 0. Rys. 3.3. Wartość funkcji celu w funkcji zmiennych (G10,G17) dla siatki o rozmiarze 6×6 węzłów (40 elementów). 0.03 0.025. f(G). 0.02 0.015 0.01 0.005 0 10 8. 10 6. 8. G33. 6. 4. 4. 2. 2 0. G10. 0. Rys. 3.4. Wartość funkcji celu w funkcji zmiennych (G10,G33) dla siatki o rozmiarze 6×6 węzłów (40 elementów). 0.1. 0.08. f(G). 0.06. 0.04. 0.02. 0 10 8. 10 6. 8. G. 28. 4. 4. 2. 2 0. G. 6. 12. 0. Rys. 3.5. Wartość funkcji celu w funkcji zmiennych (G12,G28) dla siatki o rozmiarze 6×6 węzłów (40 elementów). 20.

(21) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. 0.07 0.06. f(G). 0.05 0.04 0.03 0.02 0.01 0 10 8. 10 6. 8 6. 4. G. 4. 2. 32. 2 0. 0. G8. Rys. 3.6. Wartość funkcji celu w funkcji zmiennych (G8,G32) dla siatki o rozmiarze 6×6 węzłów (40 elementów). 0.014 0.012. f(G). 0.01 0.008 0.006 0.004 0.002 0 10 8. 10 6. 8. G29. 4. 4. 2. 2 0. G14. 6. 0. Rys. 3.7. Wartość funkcji celu w funkcji zmiennych (G14,G29) dla siatki o rozmiarze 6×6 węzłów (40 elementów). 1.4 1.2. f(G). 1 0.8 0.6 0.4 0.2 0 10 8. 10 6. G39. 8 4. 4. 2. 2 0. G26. 6. 0. Rys. 3.8. Wartość funkcji celu w funkcji zmiennych (G26,G39) dla siatki o rozmiarze 6×6 węzłów (40 elementów). 21.

(22) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. 0.8 0.7 0.6. f(G). 0.5 0.4 0.3 0.2 0.1 0 10 8. 10 6. G36. 8 6. 4. 4. 2. 2 0. G9. 0. Rys. 3.9. Wartość funkcji celu w funkcji zmiennych (G9,G36) dla siatki o rozmiarze 6×6 węzłów (40 elementów). 1.4 1.2. f(G). 1 0.8 0.6 0.4 0.2 0 10 8. 10 6. G35. 8 4. 4. 2. 2 0. G21. 6. 0. Rys. 3.10. Wartość funkcji celu w funkcji zmiennych (G21,G35) dla siatki o rozmiarze 6×6 węzłów (40 elementów). Na rys. 3.11 przedstawiono jednowymiarowe przekroje funkcji celu dla przypadku siatki o rozmiarze 10×10 węzłów, tj. wykresy funkcji: wy1. h(l) = f((1 – l)Gorg – lGtest). (3.2). gdzie Gorg jest siatką oryginalną, zaś Gtest jest rozwiązaniem uzyskanym za pomocą jednej z metod optymalizacyjnych opisanych w dalszej części rozdziału. Dla l = 0 wyrażenie (1 – l)Gorg – lGtest opisuje siatkę oryginalną i funkcja przyjmuje minimum globalne o wartości 0. Dla l = 1 wyrażenie (1 – l)Gorg – lGtest opisuje siatkę uzyskaną w wyniku zastosowania procedury optymalizacyjnej. Jak widać z załączonych wykresów, dla l = 1 obserwujemy minimum lokalne funkcji h, co oznacza, że procedura optymalizacyjna prawdopodobnie utknęła w minimum lokalnym funkcji celu. 22.

(23) Numeryczne algorytmy tomografii rezystancji siatek rezystorów W pierwszym przypadku przedstawionym na rys. 3.11 wartość funkcji h dla siatki testowej wynosi h(1) = 1,873⋅10–8, przy maksymalnym błędzie rekonstrukcji wynoszącym 0,63698. W drugim przypadku h(1) = 5,173⋅10–9, a błąd rekonstrukcji 0,36087. W pozostałych przypadkach h(1) = 3,244⋅10–7 przy błędzie 0,74729 i h(1) = 8,113⋅10–8 przy błędzie 0,61058. Oznacza to, że funkcja celu może przyjmować bardzo niewielkie wartości dla siatek testowych odległych od siatki oryginalnej, co znacznie utrudnia proces optymalizacji. Przedstawione powyżej wyniki pokazują również, że funkcja celu może posiadać wiele minimów lokalnych. −8. −7. x 10. 1.4. 1.2. 1.2. 1. 1. 0.8. 0.8. h(λ). h(λ). 1.4. 0.6. 0.6. 0.4. 0.4. 0.2. 0.2. 0. 0. 0.2. 0.4. λ. 0.6. 0.8. 1. x 10. 0. 1.2. −7. 3.5. 0.2. 0.4. 0. 0.2. 0.4. λ. 0.6. 0.8. 1. 1.2. 0.6. 0.8. 1. 1.2. −7. x 10. x 10 4. 3. 3.5 3. 2.5. 2.5. 2. h(λ). h(λ). 0. 1.5. 2 1.5. 1. 1. 0.5 0. 0.5 0. 0.2. 0.4. λ. 0.6. 0.8. 1. 0. 1.2. λ. Rys. 3.11. Przekroje jednowymiarowe funkcji celu. 3.1.3. Algorytmy metaheurystyczne Poniżej pokrótce omówiono algorytmy metaheurystyczne, które zostały zaimplementowane w procedurach rekonstrukcji konduktancji elementów siatki. Algorytmy metaheurystyczne są grupą algorytmów poszukiwania rozwiązania problemów obliczeniowych w sposób pseudolosowy. Heurystyka (z greckiego heurisko – znajduję) to zbiór specyficznych reguł, które mogą pomóc w zaproponowaniu jak najlepszego rozwiązania dla danego problemu lub grupy problemów [49]. Metaheurystyka to reguły dotyczące sposobu generowania reguł – takich, które mogą najbardziej pomóc w znalezieniu dobrego rozwiązania. Metaheurystykami nazywa się. 23.

(24) Numeryczne algorytmy tomografii rezystancji siatek rezystorów więc zasady generowania reguł dla konkretnych algorytmów heurystycznych. W ramach jednej metaheurystyki można zwykle zaproponować przynajmniej kilka algorytmów heurystycznych, będących wariantami pewnego ogólnego podejścia. Jak już wspomniano, algorytmy metaheurystyczne często inspirowane są zjawiskami fizycznymi lub przyrodniczymi. Spośród najczęściej implementowanych, warto wymienić symulowane wyżarzanie, algorytmy genetyczne i ewolucyjne, systemy mrówkowe, rój cząsteczek czy sztuczne systemy immunologiczne. W ostatnim czasie zaproponowano również algorytmy inspirowane zachowaniem robaczków świętojanskich (ang. firefly algorithm) oraz zwyczajami lęgowymi kukułek (ang. cuckoo search) [25]. Na potrzeby realizacji algorytmów rekonstrukcji wybrano trzy metody z grupy metod losowych i pseudolosowych. Zaimplementowano metodę Monte-Carlo, w której parametry nowej próbki (siatki testowej) nie zależą od parametrów próbki poprzedniej. Następnie zastosowano algorytm symulowanego wyżarzania, w którym nowa próbka jest generowana w oparciu o najlepsze dotychczasowe rozwiązanie. Ponadto sprawdzono działanie algorytmu genetycznego w podstawowej wersji, w którym w każdym kroku algorytmu analizowany jest cały zbiór (pokolenie) siatek testowych i na jego podstawie generowane jest kolejne pokolenie. Metoda Monte Carlo Algorytm ten został zaproponowany w [32] przez Nicholasa Metropolisa oraz Stanisława Ulama. Prowadzili oni badania w laboratorium w Los Alamos, poszukując metod znajdowania rozwiązań problemów matematycznych, których rozwiązanie analityczne było niemożliwe [31]. Zaproponowana przez nich metoda była w zasadzie w pełni losowa, ale dała początek całej grupie metod pseudolosowych, w których przestrzeń poszukiwań rozwiązania postawionego problemu jest przeszukiwana w uporządkowany i zalgorytmizowany sposób. Istotą metody Monte Carlo jest ocena losowo wygenerowanego rozwiązania problemu (obliczenie wartości funkcji celu) i porównanie tej wartości z wartością funkcji celu dla najlepszego dotychczas znalezionego rozwiązania. Jeżeli nowe rozwiązanie jest lepsze, to poprzednie najlepsze rozwiązanie jest usuwane i nowe jest zapisywane na jego miejsce. W metodzie tej nowe rozwiązanie jest generowane w następnym kroku algorytmu w sposób całkowicie losowy, nie biorąc pod uwagę wyników obliczeń przeprowadzonych wcześniej. Algorytm działający w oparciu o metodę Monte Carlo zrealizowano w następujący sposób. Założono, że znany jest rozmiar siatki oraz przedział w jakim mieszczą się wartości konduktancji jej elementów. W każdym kroku algorytmu. 24.

(25) Numeryczne algorytmy tomografii rezystancji siatek rezystorów generowana jest siatka o zadanym rozmiarze, w której konduktancji każdego elementu przypisywana jest wartość zmiennej losowej o rozkładzie równomiernym z zadanego przedziału. Dla tak wygenerowanej siatki obliczana jest wartość funkcji celu. Jeżeli wartość ta jest mniejsza od wartości dla dotychczas najlepszego rozwiązania, nowe rozwiązanie jest akceptowane jako najlepsze. W następnym kroku, nowa próbka jest generowana całkowicie losowo. W związku ze sporym rozmiarem problemu nie należy się spodziewać zbyt dobrych wyników działania tego algorytmu. Zostaną one zbiorczo przedstawione i porównane z wynikami osiągniętymi dla innych metod w rozdziale 4. Symulowane wyżarzanie Symulowane wyżarzanie [24] jest jednym z pierwszych algorytmów metaheurystycznych, który znalazł zastosowanie w wielu dziedzinach optymalizacji. Metoda ta polega na poszukiwaniu nowego rozwiązania w pobliżu aktualnego rozwiązania. Podstawową zaletą algorytmu symulowanego wyżarzania jest zdolność unikania utknięcia w lokalnym minimum funkcji. Prawdopodobieństwo akceptacji kolejnej próbki jest uzależnione od wartości funkcji celu próbki poprzedniej, przy czym z niezerowym prawdopodobieństwem akceptowane są rozwiązania gorsze. Wartość funkcji celu dla badanej próbki nazywana jest tutaj energią, a warunek akceptacji nowego rozwiązania jest następujący [25]: wy2. pA = e. − ΔE T. >r. (3.3). przy czym pA oznacza prawdopodobieństwo akceptacji, ΔE – różnicę energii analizowanej i poprzedniej próbki, T – temperaturę, a r – wartość losową z przedziału (0;1), którą można nazwać progiem akceptacji. Z wzoru (3.3) wynika, że lepsze rozwiązania (o niższej energii) są zawsze akceptowane (dla ΔE < 0, pA > 1), natomiast z pewnym prawdopodobieństwem, uzależnionym od temperatury T, akceptowane są gorsze rozwiązania. Ta cecha algorytmu pozwala na opuszczenie lokalnego minimum i poszukiwanie minimum globalnego. Wartość parametru T jest stopniowo obniżana w kolejnych krokach algorytmu, aż do osiągnięcia wartości zerowej. Nawiązuje to do procesu wyżarzania w metalurgii, stąd nazwa algorytmu. Obniżanie temperatury oznacza, że na początku działania algorytmu gorsze rozwiązania są akceptowane z większym prawdopodobieństwem, a w miarę wzrostu liczby iteracji prawdopodobieństwo to maleje.. 25.

(26) Numeryczne algorytmy tomografii rezystancji siatek rezystorów W implementacji algorytmu symulowanego wyżarzania do rozwiązania problemu rekonstrukcji konduktancji, nowa próbka jest generowana na podstawie najlepszej dotychczasowej próbki w ten sposób, że wartość konduktancji każdego elementu siatki jest modyfikowana o wartość losową z przedziału (–δ; δ). Parametry T oraz δ należy dobrać indywidualnie dla każdego postawionego problemu. Wybór początkowej wartości temperatury ma kluczowe znaczenie dla skuteczności działania algorytmu. Mała wartość temperatury oznacza, że gorsze rozwiązania praktycznie w ogóle nie będą akceptowane. Natomiast zbyt duża wartość T oznacza stan wysokiej energii systemu i w efekcie praktycznie losowe działanie algorytmu. Podobnie parametr δ ma istotne znaczenie. Zbyt mała wartość tego parametru będzie oznaczała wolną zbieżność metody, a co za tym idzie potrzebę wykonania dużej liczby iteracji. Z kolei zbyt duża wartość δ może spowodować zbyt duże różnice pomiędzy kolejnymi próbkami, co pogłębia losowy charakter algorytmu. Istotnym zagadnieniem związanym z implementacją algorytmu symulowanego wyżarzania jest również kontrola procesu wyżarzania, tzn. wybór odpowiedniej charakterystyki obniżania temperatury, tak aby chłodzenie systemu następowało stopniowo i osiągnęło stan zerowej energii w punkcie minimum globalnego. W literaturze zaproponowano wiele sposobów wyznaczania współczynnika chłodzenia lub charakterystyki obniżania temperatury. Najczęściej stosowane jest rozwiązanie, w którym temperatura obniżana jest w postępie geometrycznym ze współczynnikiem α z przedziału (0;1), zgodnie z poniższą zależnością: wy3. T ( t ) = T0 ⋅ α t. (3.4). gdzie T0 oznacza początkową wartość temperatury, t = 1, 2, ..., tf , a tf – maksymalną liczbę iteracji. Zaletą takiego wyboru jest, że T → 0 gdy t → ∞, a co za tym idzie nie ma potrzeby specyfikowania maksymalnej liczby iteracji, jeśli tylko wstępnie zdefiniowano z jaką dokładnością oczekiwane jest rozwiązanie. Algorytm genetyczny Opisany powyżej algorytm symulowanego wyżarzania działa w oparciu o analizę pojedynczej próbki. Przykładem algorytmu analizującego wiele próbek w każdym kroku jest algorytm genetyczny, należący do klasy algorytmów ewolucyjnych. Jest to klasyczny algorytm, którego idea działania bazuje na założeniach teorii ewolucji Darwina w biologii [23]. Działanie algorytmu można podzielić na kilka podstawowych kroków. W pierwszym kroku algorytmu, zwanym inicjalizacją, generowana jest populacja początkowa (pierwsze pokolenie), składające się z wielu osobników. W omawianym 26.

(27) Numeryczne algorytmy tomografii rezystancji siatek rezystorów przypadku jest to skończony zbiór siatek o zadanym rozmiarze i wartościach konduktancji elementów z zadanego przedziału. Następnie osobniki pierwszego pokolenia poddawane są ocenie – każdemu osobnikowi w pokoleniu przypisywana jest wartość funkcji przystosowania, obliczana na podstawie zdefiniowanej wzorem (3.1) funkcji celu. W kolejnym kroku algorytmu dokonywana jest selekcja najlepiej przystosowanych osobników. Zostaną one wprost skopiowane do następnego pokolenia lub poddane reprodukcji. Proces reprodukcji składa się z dwóch podstawowych czynności – krzyżowania i mutacji. Gdy już wygenerowane zostanie nowe pokolenie, algorytm wraca do punktu, w którym dokonywana jest ocena przystosowania wszystkich osobników [39]. Poniżej zostaną szczegółowo omówione operacje, które są podstawą działania algorytmu genetycznego: ocena, selekcja, krzyżowanie oraz mutacja. Procesy oceny i selekcji występują w algorytmach genetycznych nierozłącznie i często są definiowane jako jeden proces – selekcji. Z tego powodu procesy te zostaną opisane razem. Aby ułatwić proces selekcji, na etapie oceny należy przypisać wszystkim osobnikom w pokoleniu wartości funkcji przystosowania w taki sposób, aby osobniki najlepsze miały przypisane wartości najwyższe. Można to osiągnąć, bazując na wartości funkcji celu. Funkcja celu dla opisywanego problemu jest tak zdefiniowana, że dla próbek najbardziej zbliżonych do oryginału przyjmuje wartości najniższe. W związku z tym funkcję przystosowania zdefiniowano następująco: wy4. f p (G ) = ( x ⋅ f max − f (G)) y. (3.5). gdzie fmax oznacza wartość funkcji celu dla najgorszego osobnika w pokoleniu, zaś x ¥ 1, y > 0 są współczynnikami rzeczywistymi, pozwalającymi na regulację przebiegu funkcji fp. Po obliczeniu wartości funkcji przystosowania dla wszystkich osobników w danym pokoleniu, dokonywana jest normalizacja tych wartości, tak aby ich suma wynosiła 1. Dzięki temu, w następnej kolejności można dokonać selekcji poprzez losowanie ruletkowe. Każdemu osobnikowi przypisywany jest wycinek koła o obwodzie 1, którego długość krawędzi odpowiada wartości znormalizowanej funkcji przystosowania. Wygenerowanie liczby losowej z przedziału (0;1) oznacza wskazanie na odpowiedni wycinek koła, a co za tym idzie konkretnego osobnika w pokoleniu. Tak realizowana procedura selekcji pozwala wybierać do reprodukcji osobniki lepiej przystosowane z dużym prawdopodobieństwem, ale również z niezerowym prawdopodobieństwem osobniki gorzej przystosowane. Wartość współczynnika x większa od 1 zapewnia niezerowe prawdopodobieństwo wylosowania najgorszego osobnika. Współczynnik y pozwala regulować kształt funkcji przystosowania. Jeżeli większość osobników ma obliczoną podobną wartość funkcji celu, ustawienie współczynnika y na wartość silnie większą od 1 pozwala na zwiększenie. 27.

(28) Numeryczne algorytmy tomografii rezystancji siatek rezystorów prawdopodobieństwa wyselekcjonowania osobników lepiej przystosowanych. Jeżeli z kolei występują bardzo duże wahania wartości funkcji celu w pokoleniu, to wartość y mniejsza od 1 pozwoli na wyrównanie szans, czyli zwiększenie prawdopodobieństwa wylosowania osobników słabiej przystosowanych. Na wyselekcjonowanej grupie osobników dokonuje się procesu reprodukcji, składającego się z krzyżowania oraz mutacji. Tak powstałe osobniki tworzą kolejne pokolenie. Krzyżowanie polega na połączeniu części genotypu jednego osobnika z pozostałą częścią genotypu innego, natomiast mutacja polega na zmianie wartości losowego genu o losową wartość. Warto zwrócić uwagę, że krzyżowaniem i mutacją rządzą procesy losowe. Dla każdego osobnika, decyzja o tym czy weźmie on udział w krzyżowaniu i mutacji jest podejmowana na podstawie wartości zmiennej losowej z przedziału (0;1). Przyjęto, że krzyżowanie następuje ze stosunkowo dużym prawdopodobieństwem (typowo 0,8–0,95), natomiast mutacja z małym prawdopodobieństwem (poniżej 0,1). W przypadku krzyżowania, osobniki są dobierane w pary w sposób losowy. Nie wykluczono wylosowania do krzyżowania dwa razy tego samego osobnika. Dzięki temu osobniki najlepiej przystosowane z dużym prawdopodobieństwem są kopiowane w całości do pokolenia potomnego. W przypadku zastosowania algorytmu genetycznego do rekonstrukcji siatek rezystancyjnych przyjęto, że genotypem jest wektor konduktancji wszystkich elementów siatki, a genem wartość konduktancji pojedynczego elementu. Krzyżowanie w tym przypadku polega na wylosowaniu węzła siatki, który będzie punktem granicznym krzyżowania. Nowa siatka jest generowana poprzez skopiowanie wszystkich elementów jednej z krzyżowanych siatek do wybranego węzła, oraz elementów drugiej siatki za wybranym węzłem. W procesie mutacji losowana jest liczba elementów jej podlegających, ich położenie, a następnie przypisywana jest im całkowicie nowa wartość konduktancji z ustalonego przedziału. Zaobserwowano, że w badanych przypadkach, gdy zrezygnowano z mutacji (prawdopodobieństwo mutacji ustalono na 0,0) to po określonej liczbie generacji (zależnie od rozmiaru siatki) cała populacja składała się z osobników o takim samym genotypie. 3.1.4. Algorytmy optymalizacyjne Opisane powyżej algorytmy metaheurystyczne często są w literaturze zaliczane do grupy szeroko rozumianych algorytmów optymalizacyjnych. W niniejszej pracy postanowiono jednak rozdzielić te dwie grupy algorytmów. Wyszczególniono grupę algorytmów pseudolosowych, mających swą genezę w obserwacji procesów fizycznych i przyrodniczych oraz grupę algorytmów, których zasada działania opiera się na wiedzy matematycznej z zakresu właściwości funkcji wielu zmiennych oraz metod. 28.

(29) Numeryczne algorytmy tomografii rezystancji siatek rezystorów poszukiwania ich minimów. Druga z wymienionych grup została w niniejszej pracy określona jako grupa algorytmów optymalizacyjnych i jest opisana poniżej. Wybrane algorytmy optymalizacyjne można podzielić na dwie grupy – takie, które używają bądź nie używają obliczonych numerycznie lub analitycznie wartości pochodnych funkcji celu do poszukiwania jej minimum. Jednym z algorytmów nie korzystających z informacji dostarczanych przez pochodne funkcji celu jest algorytm sympleksowy Neldera-Meada. Spośród algorytmów korzystających z pochodnych wybrano metodę quasi-Newtona z metodą interpolacji wielomianowej, metodę sekwencyjnego programowania kwadratowego oraz nieliniową metodę najmniejszych kwadratów. Algorytm sympleksowy Neldera-Meada Metody optymalizacyjne wykorzystujące wartości pochodnych funkcji celu mogą być bardzo efektywne, ale mogą też stawiać bardzo restrykcyjne wymagania co do własności funkcji celu. Uzyskanie analitycznie bądź numerycznie pochodnych funkcji celu często nie jest trywialne i wiąże się ze znacznym kosztem obliczeniowym. Ponadto w przypadku, gdy funkcja celu wykazuje silną nieciągłość, metody nie używające pochodnych wydają się bardziej naturalne [25]. Algorytm sympleksowy Neldera-Meada to metoda bezpośredniego poszukiwania, nie używająca pochodnych funkcji celu [26]. Sympleksem nazywana jest bryła o niezerowej objętości rozpięta na p+1 wierzchołkach w przestrzeni p-wymiarowej. Taka liczba wierzchołków pozwala na przeszukiwanie przestrzeni p-wymiarowej we wszystkich możliwych kierunkach. W pierwszej kolejności generowany jest sympleks wokół punktu początkowego x0. Istnieje wiele sposobów wyznaczania punktów sympleksu startowego. W badaniach opisanych w niniejszej pracy zastosowano następującą regułę. Jednym z wierzchołków sympleksu jest punkt x0, natomiast pozostałe p wierzchołków powstaje poprzez dodanie 5% wartości do każdej współrzędnej punktu x0. W kolejnych krokach algorytm modyfikuje sympleks wielokrotnie, na podstawie procedury opisanej poniżej [26]. Pozwala to na zmierzanie w uporządkowany sposób do minimum globalnego funkcji. W przedstawionej procedurze zostaną wyróżnione operacje odbicia, ekspansji, kontrakcji i redukcji, które wykonywane w odpowiedniej sekwencji realizują procedurę Neldera-Meada. Należy zwrócić uwagę, że nie wszystkie przedstawione poniżej kroki muszą być wykonane w każdej iteracji algorytmu. W większości kroków zdefiniowany jest warunek po osiągnięciu którego rozpoczynana jest kolejna iteracja.. 29.

(30) Numeryczne algorytmy tomografii rezystancji siatek rezystorów Algorytm Neldera-Meada: 1. Oznacz punkty sympleksu przez xi, gdzie i = 1, ..., p+1. 2. Uporządkuj punkty sympleksu względem wartości funkcji celu od najmniejszej f(x1) do największej f(xp+1). W każdym kroku danej iteracji, algorytm odrzuca obecnie najgorszy punkt xp+1 i dodaje nowy punkt do sympleksu (lub w przypadku kroku nr 7 zmienia wszystkie p punktów sympleksu poza x1). 3. Wygeneruj punkt odbity xr = 2 ⋅ xm − x p +1 , gdzie 1 p xm = ⋅ ∑ xi wy5 p i =1. jest środkiem ciężkości sympleksu. Następnie dla powstałego w ten sposób punktu odbitego oblicz wartość funkcji celu f(xr). 4. Jeżeli f ( x1 ) ≤ f ( xr ) ≤ f ( x p ) to zaakceptuj punkt xr i zakończ tę iterację (odbicie). 5. Jeżeli f ( xr ) ≤ f ( x1 ) to wygeneruj punkt ekspansji xs = xm + 2 ⋅ ( xm − x p +1 ) i oblicz wartość funkcji celu w tym punkcie f(xs). a. Jeżeli f ( xs ) < f ( xr ) to zaakceptuj punkt xs i zakończ tę iterację (ekspansja). b. W przeciwnym razie zaakceptuj punkt xr i zakończ tę iterację (odbicie). 6. Jeżeli f ( xr ) ≥ f ( x p ) to wykonaj kontrakcję pomiędzy punktami xm oraz lepszym z punktów xp+1 i xr: 1 a. Jeżeli f ( xr ) < f ( x p +1 ) to wygeneruj punkt kontrakcji xc = xm + ⋅ ( xr − xm ) 2 i oblicz f(xc). Jeżeli f ( xc ) < f ( xr ) to zaakceptuj xc i zakończ tą iterację (kontrakcja na zewnątrz). W przeciwnym razie przejdź do kroku 7. 1 b. Jeżeli f ( xr ) ≥ f ( x p +1 ) to wygeneruj punkt kontrakcji xcc = xm + ⋅ ( x p +1 − xm ) 2 i oblicz f(xcc). Jeżeli f ( xcc ) < f ( x p +1 ) to zaakceptuj xcc i zakończ tą iterację (kontrakcja do wewnątrz). W przeciwnym razie przejdź do kroku 7. 1 7. Wygeneruj p punktów vi = x1 + ⋅ ( xi − x1 ) i oblicz f(vi) dla i = 2, ..., p+1. 2 Sympleks kolejnej iteracji składa się z punktów x1, v2, ..., vp+1 (redukcja).. Na rys. 3.12 przedstawiono przykładowy zestaw wszystkich punktów, które mogą być wyznaczone podczas jednego kroku algorytmu Neldera-Meada w przestrzeni dwuwymiarowej. Punkty oryginalnego sympleksu połączone są linią ciągłą. Kolejne iteracje są wykonywane aż do momentu spełnienia warunków kryterium zatrzymania algorytmu.. 30.

(31) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. Rys. 3.12. Konstrukcja punktów sympleksu Neldera-Meada w przestrzeni dwuwymiarowej. Poniżej zostały pokrótce opisane kolejne algorytmy optymalizacyjne wykorzystane do rozwiązania postawionego w pracy problemu. Algorytmy te są rozbudowane i skomplikowane, dlatego zdecydowano się na opis skrócony. Szczegółowy opis metody quasi-Newtona można znaleźć w pracach [4, 11]. Algorytm sekwencyjnego programowania kwadratowego opisano między innymi w pracach [16, 20, 38]. Dokładny opis nieliniowej metody najmniejszych kwadratów znajduje się w pracach [6, 7].. Metoda quasi-Newtona. Opisana poniżej metoda quasi-Newtona z procedurą interpolacji wielomianowej zaliczana jest do grupy metod korzystających z wartości pochodnych funkcji celu. Metody te zwane są również metodami gradientowymi. W metodzie tej wartości elementów hesjanu (macierzy drugich pochodnych) nie są obliczane bezpośrednio. Przybliżenia hesjanu są obliczane na podstawie kolejnych wektorów gradientu (lub przybliżeń wektora gradientu). Metody quasi-Newtona są uogólnieniem metody siecznych poszukiwania zer pierwszej pochodnej dla problemów wielowymiarowych. Dla problemów wielowymiarowych równanie siecznych nie definiuje rozwiązania w sposób jednoznaczny, a różne wersje metody quasi-Newtona różnią się sposobem nakładania dodatkowych ograniczeń na rozwiązanie. W metodzie quasi-Newtona, w każdej iteracji dostarczana jest informacja o krzywiźnie funkcji celu, pozwalająca sformułować model problemu w postaci: wy19. 1 g (d ) = f ( xk + d ) ≈ d T Hd + cT d + b 2. 31. (3.6).

(32) Numeryczne algorytmy tomografii rezystancji siatek rezystorów gdzie H oznacza hesjan funkcji f, c jest gradientem f, a b jest wartością f w punkcie xk. Optymalnym rozwiązaniem tak postawionego problemu jest punkt, w którym pochodne cząstkowe funkcji celu są równe zeru, np.: wy20. g ' (d *) = Hd * +c = 0. (3.7). Kierunek optymalnego rozwiązania, oznaczony jako d*, można zatem opisać równaniem: wy21. d * = − H −1c. (3.8). Metoda Newtona (w przeciwieństwie do metod quasi-Newtona) oblicza hesjan wprost i na tej podstawie porusza się w kierunku spadku funkcji celu. Obliczanie numeryczne hesjanu jest bardzo złożone obliczeniowo i czasochłonne. Metody quasiNewtona pozwalają unikać tej niedogodności, gdyż na podstawie obserwacji zachowania funkcji f(x) oraz jej gradientu ∇f(x) generują informacje o krzywiźnie funkcji pozwalające przybliżać wartości hesjanu przy użyciu odpowiednich technik aproksymacyjnych. Opracowano wiele metod aproksymacji wartości hesjanu. Wykorzystana w niniejszej pracy metoda BFGS (Broyden-Fletcher-Goldfarb-Shanno) [4, 18] jest uważana za najbardziej skuteczną i uniwersalną. W metodzie tej kolejne aproksymacje hesjanu wyznaczane są na podstawie zależności: T. wy22. H k +1 = H k +. T. T. qk qk H s sH − k Tk k k T qk sk sk H k sk. (3.9). gdzie sk = xk +1 − xk , qk = ∇f ( xk +1 ) − ∇f ( xk ) . Jako punkt startowy H0 może być wybrana dowolna macierz symetryczna, dodatnio określona. Przykładowo, może to być macierz jednostkowa I. W celu uniknięcia konieczności obliczania macierzy odwrotnej Hk–1 oblicza się bezpośrednio kolejne przybliżenia macierzy Hk–1 stosując wzór ShermanaMorrisona. Wartości gradientu funkcji celu mogą być obliczone analitycznie lub przy użyciu pochodnych cząstkowych obliczonych numerycznie, za pomocą metody różnic skończonych. To wymaga zaburzenia wartości każdej zmiennej funkcji celu i obliczenia ilorazu zmiany wartości funkcji celu i wartości zaburzenia. W każdej iteracji, po wyznaczeniu macierzy Hk–1, wykonywana jest procedura interpolacji wielomianowej w kierunku: wy24. −1. d k = − H k ⋅ ∇f ( xk ). (3.10). Procedura ta wyznacza kolejny punkt xk+1 na prostej zawierającej obecny punkt xk, równoległej do kierunku przeszukiwania dk:. 32.

(33) Numeryczne algorytmy tomografii rezystancji siatek rezystorów. wy25. xk +1 = xk + αd k. (3.11). gdzie współczynnik α jest wielkością skalarną, oznaczającą długość kroku poszukiwania. Procedura interpolacji wielomianowej znajduje punkt o niższej wartości funkcji celu, znajdujący się na prostej opisanej wzorem (3.11). Funkcja celu jest interpolowana za pomocą wielomianu. Następnie wyznaczane jest minimum wielomianu interpolującego funkcję celu w wybranym obszarze. Procedura ta ma trzy podstawowe kroki: 1. Wyznaczenie odcinka, z którego będzie pochodził punkt kolejnej iteracji głównego algorytmu, tzn. wyznaczenie przedziału do jakiego będzie należał współczynnik α. 2. Wyznaczenie wielomianu interpolującego funkcję celu w tym przedziale. 3. Obliczenie minimum wielomianu interpolującego. Współczynnik α jest dobierany w ten sposób, aby były spełnione tzw. warunki Wolfe’a: wy26. f ( xk + αd k ) ≤ f ( xk ) + c1α∇f ( xk )T d k. (3.12). wy27. ∇f ( xk + αd k )T d k ≥ c2∇f ( xk )T d k. (3.13). gdzie c1 oraz c2 są stałymi, takimi że: 0 < c1 < c2 < 1. Warunek (3.12) zapewnia wystarczające obniżenie wartości funkcji celu. Warunek (3.13) zdefiniowano w taki sposób, aby długość kroku poszukiwania nie była zbyt mała. Punkty, które spełniają oba powyższe warunki określa się jako punkty dopuszczalne. Sekwencyjne programowanie kwadratowe. Kolejnym zastosowanym algorytmem optymalizacyjnym jest metoda sekwencyjnego programowania kwadratowego, należąca do klasy metod optymalizacji nieliniowej z ograniczeniami [16]. Stosowanie metod optymalizacji z ograniczeniami do rozwiązania omawianego w pracy problemu jest uzasadnione. Wynika to z założenia, że na początku procedury rekonstrukcji znany jest zakres, w którym mieszczą się wartości konduktancji wszystkich elementów siatki. W tej klasie algorytmów podstawowym celem jest przekształcenie problemu w łatwiejszy do rozwiązania „podproblem” bez ograniczeń, który może posłużyć za bazę procesu iteracyjnego. Przekształcenia takiego dokonuje się poprzez użycie funkcji wykluczającej punkty będące poza brzegiem ograniczenia. W ten sposób problem. 33.

(34) Numeryczne algorytmy tomografii rezystancji siatek rezystorów optymalizacji z ograniczeniami jest rozwiązywany przy użyciu sekwencji sparametryzowanych optymalizacji bez ograniczeń, które w granicy zmierzają do rozwiązania problemu z ograniczeniami. W problemie poszukiwania minimum funkcji celu f(x) metodą optymalizacji z ograniczeniami, ograniczenia mają postać: wy28. gi(x) = 0 ,. dla i = 1, ..., me. (3.14). wy29. gi(x) ≤ 0 ,. dla i = me+1, ..., m. (3.15). Ograniczenia (3.14) i (3.15) są nazywane odpowiednio ograniczeniami typu równościowego i nierównościowego. W celu poprawy efektywności algorytmów optymalizacji z ograniczeniami, stosuje się metody poszukujące rozwiązania równań Karush’a-Kuhn’a-Tucker’a (KKT). Jeżeli problem jest tzw. problemem programowania wypukłego, tzn. f(x) oraz gi(x) są funkcjami wypukłymi, wtedy równania KKT stanowią warunek konieczny i wystarczający istnienia globalnego rozwiązania problemu minimalizacyjnego. W odniesieniu do problemu poszukiwania minimum funkcji celu f(x) równania te mogą być sformułowane następująco: m. wy30. ∇f ( x*) + ∑ λi ⋅ ∇gi ( x*) = 0. (3.16). i =1. wy31. λi ⋅ gi ( x*) = 0 , dla i = 1, ..., me λi ≥ 0 , dla i = me+1, ..., m. wy32. (3.17) (3.18). gdzie x* oznacza optymalne rozwiązanie problemu. Równanie (3.16) opisuje redukcję gradientów w punkcie optymalnego rozwiązania dla funkcji f(x) oraz aktywnych ograniczeń. Dla tych gradientów, współczynniki Lagrange’a (λi, i = 1, ..., m) muszą być tak dobrane, aby zbilansować różnicę wartości funkcji celu oraz gradientów ograniczenia. W operacji redukcji mogą brać udział tylko aktywne ograniczenia, więc dla tych nieaktywnych współczynniki Lagrange’a ustawione są na 0. Warunek ten wyrażony jest zależnościami (3.17) i (3.18). Ograniczenie nierównościowe gi(x) ≤ 0 nazywane jest aktywnym w punkcie dopuszczalnym xˆ , jeżeli gi ( xˆ ) = 0 . Gdy spełniony jest warunek gi ( xˆ ) < 0 , ograniczenie nazywamy nieaktywnym. Ograniczenia w postaci równań zawsze występują w zbiorze ograniczeń aktywnych [19]. Rozwiązanie równań KKT stanowi podstawę wielu algorytmów optymalizacji nieliniowej. W metodach quasi-Newtona z ograniczeniami, wartości przybliżeń hesjanu spełniające równania KKT obliczane są przy użyciu odpowiednich procedur aktualizacyjnych, np. formuły BFGS. Metody te nazywane są metodami sekwencyjnego 34.

(35) Numeryczne algorytmy tomografii rezystancji siatek rezystorów programowania kwadratowego (ang. Sequential Quadratic Programming – SQP), gdyż w każdej iteracji głównego algorytmu rozwiązywany jest podproblem programowania kwadratowego. Główną ideą opisywanej metody dla problemów minimalizacyjnych jest sformułowanie podproblemu programowania kwadratowego w oparciu o aproksymację kwadratową funkcji Lagrange’a: m. wy33. L( x, λ ) = f ( x) + ∑ λi ⋅ g i ( x). (3.19). i =1. Warunki KKT przyjmują wtedy postać: wy34. 1 minn d T H k d + ∇f ( xk )T d d∈ℜ 2. (3.20). wy35. ∇gi ( xk )T d + gi ( xk ) = 0 , dla i = 1, ..., me. (3.21). wy36. ∇gi ( xk )T d + gi ( xk ) ≤ 0. (3.22). , dla i = me+1, ..., m. Rozwiązanie takiego podproblemu jest wykorzystywane do sformułowania kolejnej iteracji rozwiązania. W tym przypadku macierz Hk jest dodatnio określoną aproksymacją hesjanu funkcji Lagrange’a. Problem minimalizacji zdefiniowany jako problem z ograniczeniami nieliniowymi często może być rozwiązany w mniejszej liczbie iteracji niż w przypadku, gdy zostanie zdefiniowany jako problem bez ograniczeń. Wynika to z faktu, iż dzięki ograniczeniu przestrzeni poszukiwań algorytm może w sposób bardziej efektywny określać kierunek i długość kroku przeszukiwania. Implementacja algorytmu sekwencyjnego programowania kwadratowego zawiera kilka podstawowych kroków. W pierwszej kolejności wykonywane jest obliczenie przybliżenia hesjanu. Aproksymacja hesjanu funkcji Lagrange’a Hk jest obliczana metodą BFGS. Aby hesjan był dodatnio określony należy zapewnić, że wartość qkTsk jest dodatnia dla każdego przybliżenia oraz zainicjować aproksymację macierzy H macierzą dodatnio określoną. W każdej iteracji głównego algorytmu sekwencyjnego programowania kwadratowego rozwiązywany jest podproblem programowania kwadratowego. Procedura wyznaczenia rozwiązania podproblemu programowania kwadratowego składa się z dwóch faz. Pierwsza faza obejmuje wyznaczenie punktu dopuszczalnego. W drugiej fazie generowany jest ciąg punktów dopuszczalnych zbieżnych do rozwiązania. W tym celu stosowana jest metoda zbioru aktywnych ograniczeń (ang. active set) [15, 16]. W pierwszej części algorytmu wyznacza się punkt startowy, rozwiązując problem programowania liniowego. Następnie może być zrealizowana główna część 35.