Index of /rozprawy2/10700

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. Kształtowanie rozkładu pola magnetycznego z wykorzystaniem zaawansowanych metod obliczeniowych. mgr inż. Bartłomiej Garda. Promotor: dr hab. inż. Zbigniew Galias, prof. AGH. Kraków, 2013.

(2) Spis treści. Spis treści. i. Wykaz ważniejszych skrótów i oznaczeń. iv. 1 Wstęp 1.1 Cel i teza pracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Zawartość pracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Przedstawienie problemu 2.1 Wprowadzenie . . . . . . . . . . . . 2.2 Synteza pola magnetycznego . . . . . 2.3 Pole magnetyczne cewki . . . . . . . 2.4 Przypadek liniowy . . . . . . . . . . 2.5 Błąd dyskretyzacji . . . . . . . . . . 2.5.1 Cewka nieskończenie cienka . 2.5.2 Cewka o skończonej grubości 2.6 Definicja funkcji celu . . . . . . . . . 2.7 Ocenianie wyników . . . . . . . . . .. 1 2 3. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. 4 4 5 6 9 10 10 11 12 14. 3 Metody optymalizacji 3.1 Podstawowe definicje . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Optymalizacja z ograniczeniami . . . . . . . . . . . . . . . . . . 3.3 Klasyfikacja metod optymalizacji . . . . . . . . . . . . . . . . . 3.4 Zagadnienie najmniejszych kwadratów . . . . . . . . . . . . . . 3.4.1 Zagadnienie najmniejszych kwadratów bez ograniczeń . 3.4.2 Zagadnienie najmniejszych kwadratów z ograniczeniami 3.5 Metoda regularyzacji . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Optymalizacja nieliniowa . . . . . . . . . . . . . . . . . . . . . . 3.6.1 Metoda sympleksowa Neldera-Meada . . . . . . . . . . . 3.6.2 Gradientowe metody optymalizacji bez ograniczeń . . . 3.6.3 Metody gradientowe dla optymalizacji z ograniczeniami. . . . . . . . . . . .. . . . . . . . . . . .. 16 16 17 19 21 21 24 30 32 32 33 38. i. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . ..

(3) ii. SPIS TREŚCI. 3.7. 3.6.4 Nieliniowe zagadnienie najmniejszych kwadratów Algorytmy niedeterministyczne . . . . . . . . . . . . . . 3.7.1 Algorytm symulowanego wyżarzania . . . . . . . 3.7.2 Algorytmy ewolucyjne . . . . . . . . . . . . . . .. . . . .. . . . .. . . . .. . . . .. 4 Optymalizacja rozkładu prądu w cewce 4.1 Opis rozważanych problemów . . . . . . . . . . . . . . . . . . . 4.2 Optymalizacja pola magnetycznego bez ograniczeń . . . . . . . 4.2.1 Wybór metody rozwiązania układu równań normalnych 4.2.2 Pole magnetyczne kształtowane na osi . . . . . . . . . . 4.2.3 Pole magnetyczne kształtowane na powierzchni kuli . . 4.2.4 Podsumowanie problemu optymalizacji bez ograniczeń . 4.3 Regularyzacja Tichonowa . . . . . . . . . . . . . . . . . . . . . 4.4 Optymalizacja z ograniczeniami . . . . . . . . . . . . . . . . . . 4.4.1 Nieujemna metoda najmniejszych kwadratów . . . . . . 4.4.2 Zagadnienie najmniejszych kwadratów z ograniczeniami 4.5 Porównanie wyników . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Inne metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.1 Metoda quasi-Newtona . . . . . . . . . . . . . . . . . . 4.6.2 Metoda sympleksowa Neldera–Meada . . . . . . . . . . 4.6.3 Symulowane wyżarzanie . . . . . . . . . . . . . . . . . . 4.6.4 Algorytm genetyczny . . . . . . . . . . . . . . . . . . . . 4.6.5 Podsumowanie wyników . . . . . . . . . . . . . . . . . . 4.7 Synteza niejednorodnego pola magnetycznego . . . . . . . . . . 5 Optymalizacja kształtu cewki 5.1 Przedstawienie problemu . . . . . . . . . . . . . . . . . . . . 5.2 Funkcja celu dla problemu optymalizacji kształtu . . . . . . 5.3 Wyniki obliczeń . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Pole magnetyczne kształtowane na osi . . . . . . . . 5.3.2 Podsumowanie wyników . . . . . . . . . . . . . . . . 5.3.3 Optymalizacja kształtu cewki z gęstością prądu . . . 5.3.4 Inne metody . . . . . . . . . . . . . . . . . . . . . . 5.3.5 Pole magnetyczne kształtowane na powierzchni kuli 5.3.6 Podsumowanie wyników . . . . . . . . . . . . . . . . 5.4 Synteza niejednorodnego pola magnetycznego . . . . . . . .. . . . . . . . . . .. . . . . . . . . . .. . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . .. . . . .. 40 41 41 42. . . . . . . . . . . . . . . . . . .. 48 48 50 50 52 55 59 62 67 70 76 81 82 85 87 87 88 89 90. . . . . . . . . . .. 95 95 97 98 98 101 106 109 112 114 115. 6 Analiza rzeczywistej cewki. 121. 7 Podsumowanie. 125. Bibliografia. 129.

(4) SPIS TREŚCI. A Pole magnetyczne zwoju kołowego. iii. 135.

(5) Wykaz ważniejszych skrótów i oznaczeń. Γ. –. Ω m n I, x j A B, Br , Bz breq , b f g K(·), E(·) K(·, ·, ·, ·) L λ, µ. – – – – – – – – – – – – – –. µ0 = 4π10−7 ζ(z) rc rcmin rcmax zcmin zcmax rreq zreqmin. – – – – – – – – –. zreqmax. –. obszar w którym generowane jest pole magnetyczne o zadanym rozkładzie obszar obejmujący przewody cewki liczba punktów w obszarze Γ liczba punktów w obszarze Ω wektor natężeń prądu elektrycznego w uzwojeniach gęstości prądu w uzwojeniu cewki potencjał wektorowy pola magnetycznego wektor pola magnetycznego oraz jego składowe składowa w kierunku z zadanego pola magnetycznego funkcja celu funkcja ograniczeń całki eliptyczne, pierwszego i drugiego rodzaju jądro równania całkowego Fredholma funkcja Lagrange’a współczynniki Lagrange’a odpowiadające ograniczeniom równościowym i nierównościowym przenikalność magnetyczna próżni funkcja opisująca wewnętrzny kształt cewki promień cewki w przypadku cewki jednowymiarowej wewnętrzny promień cewki w przypadku dwuwymiarowym zewnętrzny promień cewki w przypadku dwuwymiarowym współrzędna z początku cewki współrzędna z końca cewki promień sfery z zadanym polem magnetycznym współrzędna z początku obszaru z zadanym polem magnetycznym współrzędna z początku obszaru z zadanym polem magnetycznym. iv.

(6) Podziękowania. Serdecznie chciałbym podziękować prof. dr. hab. inż Zbigniewowi Galiasowi, mojemu promotorowi za pomoc, zaangażowanie i mnóstwo cennych wskazówek, bez których z pewnością niniejsza praca nie mogłaby powstać. Składam również podziękowania prof. dr. hab. inż. Antoniemu Cieśli za wskazanie i wprowadzenie w tak ciekawy temat, jakim jest synteza pola magnetycznego.. Pracę dedykuję moim synom Piotrowi i Michałowi.. v.

(7)

(8) Rozdział 1. Wstęp. Budowa magnesów wytwarzających pole magnetyczne o zadanym kształcie należy do zagadnień syntezy pola magnetycznego. Synteza pola jest przedmiotem intensywnych badań od wielu lat. Już w latach 50. pojawiły się pierwsze prace dotyczące kształtowania pola magnetycznego [6][33]. W latach 60. opublikowano ważne prace dotyczące konstrukcji cewek rezystancyjnych i nadprzewodzących [5][31][17][59]. Zastosowano aparat optymalizacji analitycznej wykorzystującej pewne uproszczone modele. Ponieważ, w ogólności synteza pola magnetycznego jest zagadnieniem źle postawionym w sensie Hadamarda [39], ważną rolę w badaniach nad syntezą pola magnetycznego odegrało wykorzystanie metod regularyzacji, szczególnie regularyzacji Tichonowa [1][23][75][74][67]. Ze względu na fakt liniowej zależności pomiędzy prądem a wytworzonym przez niego polem magnetycznym pojawiło się wiele rozwiązań wykorzystujących programowanie liniowe oraz metodę najmniejszych kwadratów [41][57][62]. W wielu pracach dotyczących kształtowania pola magnetycznego wykorzystuje się metody niedeterministyczne. Wykorzystane są np. metody Tabu-Search [56] oraz algorytmy ewolucyjne i genetyczne [20][50][71]. Jednym z najczęściej pojawiających się problemów syntezy pola magnetycznego jest wytworzenie w pewnym zadanym obszarze pola jednorodnego. Problem ten znajduje zastosowanie przy budowie urządzeń na potrzeby diagnostyki medycznej takich jak rezonans magnetyczny (Magnetic Resonance Imaging, MRI) [46][69], gdzie silne jednorodne pole musi być wzbudzone w założonym obszarze. Autor w swojej dotychczasowej pracy zetknął się z problemem kształtowania rozkładu pola w zastosowaniu magnesów do budowy filtrów magnetycznych [14][66]. W tym przypadku wytworzenie jednorodnego pola magnetycznego stwarza korzystne warunki do ekstrakcji cząstek z filtrowanej zawiesiny[15]. Zarówno w urządzeniach MRI jak i w filtrach magnetycznych wymagane jest wzbudzanie silnych pól magnetycznych, w obszarach o znacznej objętości. Najczęściej wykorzystuje się do tego celu uzwojenia nadprzewodnikowe. 1.

(9) 2. 1.1. ROZDZIAŁ 1. WSTĘP. Cel i teza pracy. W początkowym okresie swojej pracy badawczej autor zetknął się z problemem rozpoznania kształtu i parametrów wewnętrznych magnesu nadprzewodnikowego znajdującego się Laboratorium Kriotechniki Katedry Elektrotechniki i Elektorenergetyki. Na podstawie pomiarów pośrednich (w tym pomiaru pola magnetycznego wewnątrz cewki) dokonano zgrubnej analizy kształtu cewki wraz z dozwojeniami. Następnie autor dokonał próby optymalizacji kształtu pola magnetycznego z wykorzystaniem algorytmów genetycznych oraz Metody Elementów Skończonych (MES) [28][30]. Podejście takie nie dało zadowalających rezultatów szczególnie dla większych wymiarów problemu. Dlatego w niniejszej dysertacji autor postawił następującą tezę: Wykorzystanie odpowiednio dobranych metod obliczeniowych, umożliwia dokładniejsze i bardziej efektywne rozwiązanie zagadnienia syntezy pola magnetycznego, niż metody dotychczas stosowane. Teza ta wynikała z pewnych założeń dotyczących analizowanego problemu. Rezygnując w pewnym sensie z ogólności rozwiązania, można konstruować cewki wytwarzające zadane pole magnetyczne w sposób bardziej dokładny. Strata ogólności problemu polega na rezygnacji z metody elementów skończonych i wykorzystaniu analitycznych zależności pomiędzy rozważanymi wielkościami. Rezygnacja z metody elementów skończonych była możliwa, gdyż rozważano problem wytworzenia pola magnetycznego w środowisku magnetycznie liniowym (powietrze). W celu potwierdzenia postawionej tezy wyznaczono następujący cel pracy: Celem pracy jest porównanie różnych metod obliczeniowych do rozwiązania problemu syntezy pola magnetycznego oraz wskazanie najbardziej efektywnej metody. Metody będą porównane pod względem jakości rozwiązania oraz złożoności obliczeniowej. Analizowane będą liniowe oraz nieliniowe problemy syntezy pola, tzn. problem doboru wartości prądów przy zadanym kształcie cewki jak i problem doboru kształtu cewki przy stałej gęstości prądu. Pole magnetyczne możemy kształtować za pomocą zmiany wartości prądów w przewodach, których położenie jest ustalone lub poprzez zmianę położenia przewodów i ewentualnie również wartości prądów w tych przewodach. Ponieważ w środowiskach magnetycznie liniowych związek pomiędzy wartością prądu a polem przez niego wytworzonym jest liniowy, to pierwszy przypadek kształtowania pola (dobór wartości prądów przy ustalonym położeniu) jest problemem liniowym. Ponieważ położenie przewodu w sposób nieliniowy wpływa na wartość pola magnetycznego, to drugi przypadek kształtowania pola nazywamy.

(10) 1.2. ZAWARTOŚĆ PRACY. 3. nieliniowym.. 1.2. Zawartość pracy. Praca składa się z wykazu najważniejszych oznaczeń, siedmiu rozdziałów, dodatku oraz spisu literatury. Pierwszy rozdział stanowi wstęp oraz przegląd literatury. Postawiono w nim tezę pracy oraz nakreślono jej cel. Streszczono również zawartość pracy. W rozdziale drugim zdefiniowano problem syntezy pola magnetycznego. Dokonano analizy błędów dyskretyzacji. Zdefiniowano funkcję celu dla problemu optymalizacyjnego, oraz zdefiniowano wielkości używane do oceny wyników. Rozdział trzeci przedstawia wybrane metody optymalizacji funkcji wielu zmiennych. W rozdziale tym dokonano przeglądu metod optymalizacji numerycznej z ograniczeniami oraz bez ograniczeń. Opisano metody wykorzystane w niniejszej pracy. W rozdziale czwartym analizowany jest problem syntezy w przypadku liniowym, w którym pole magnetyczne kształtowane jest za pomocą zmiany rozkładu prądów w cewce. Dokonano analizy wielu metod dla problemu z ograniczeniami i bez ograniczeń, porównano wyniki i wskazano najlepszą metodę. W rozdziale piątym przedstawiono analizę problemu syntezy pola poprzez zmianę kształtu cewki przy ustalonej gęstości prądu. Do rozwiązania problemu wykorzystano techniki optymalizacji nieliniowych funkcji wielu zmiennych. Dokonano porównania wyników uzyskanych różnymi metodami oraz przedstawiono wnioski. W rozdziale szóstym przedstawiono analizę rzeczywistej cewki, wykonanej wg projektu autora. Dokonano porównania z cewką zaprojektowaną przy użyciu metod analizowanych w niniejszej pracy. Podsumowanie dysertacji zawarto w rozdziale siódmym. Zostały tam przedstawione wnioski wyciągnięte na podstawie przeprowadzonych badań. Zawarto podsumowanie zrealizowanych celów oraz ustosunkowanie się autora do postawionej tezy pracy. W dodatku przedstawiono wyprowadzenie analitycznych zależności pomiędzy pojedynczym zwojem prądowym a polem magnetycznym w ustalonym punkcie przestrzeni..

(11) Rozdział 2. Przedstawienie problemu. W rozdziale tym przedstawiony zostanie problem kształtowania pola w cewce powietrznej. Problem ten zostanie sformułowany jako problem optymalizacji.. 2.1. Wprowadzenie. Problemy współczesnej techniki często dzieli się na problemy analizy oraz problemy syntezy. Analiza polega na wyznaczeniu odpowiedzi układu na zadane wymuszenie. Problemy takie opisane są za pomocą pewnych matematycznych zależności, które w odpowiedni sposób przybliżają (modelują) dane zagadnienie. Problemy te opisane są zwykle za pomocą równań algebraicznych i równań różniczkowych zwyczajnych lub cząstkowych. Równania te konstruuje się wraz z pewnymi warunkami początkowymi lub granicznymi. Dzięki rozwojowi technik komputerowych oraz metod numerycznych rozwiązanie problemów analizy zwykle nie sprawia większych problemów. Niniejsza praca zajmuje się problemami odwrotnymi do problemów analizy. Problemy takie nazywamy problemami syntezy. W zagadnieniach tych poszukiwane jest wymuszenie, które działając na dany układ spowoduje odpowiednią (zadaną z góry) odpowiedź. Ogólnie związek pomiędzy wymuszeniem a odpowiedzią można opisać równaniem: F {x} = y. (2.1). Równanie (2.1) z operatorem F przedstawia problem analizy, gdy poszukiwana jest odpowiedź układu y, zaś znajdowanie wymuszenia x to problem syntezy. Operator F może być operatorem liniowym lub nieliniowym, stacjonarnym lub niestacjonarnym. Problem syntezy w ogólnym przypadku może nie posiadać rozwiązania lub może nie mieć własności jednoznaczności rozwiązań. Przykładem wziętym z życia może być problem rekonstrukcji obrazu, gdzie wejściem jest fotografowany obiekt, systemem jest układ optyczny aparatu, natomiast wyjściem 4.

(12) 2.2. SYNTEZA POLA MAGNETYCZNEGO. 5. są dane otrzymane przez przetwornik optyczny. Problemem syntezy byłoby odtworzenie obiektu na podstawie jego zdjęcia. Jest to naturalny przykład sytuacji o niejednoznacznym rozwiązaniu problemu syntezy.. Rysunek 2.1: Problem analizy i syntezy. 2.2. Synteza pola magnetycznego. W niniejszej pracy będzie analizowany problem syntezy pola magnetycznego. Polega on na znalezieniu takiego rozkładu gęstości prądów w pewnym obszarze, aby pole magnetyczne wytworzone przez ten prąd było takie, jak założono. Można tutaj ”sterować” zarówno wartością prądu jak również rozmieszczeniem przewodów w przestrzeni. W trakcie przeprowadzonych badań nad problemem syntezy pola magnetycznego dokonano pewnych uproszczeń. Analizę ograniczono do statycznego pola magnetycznego oraz założono, że pole magnetyczne jest wytworzone w materiałach o liniowych, jednorodnych i izotropowych właściwościach magnetycznych. Założono ponadto, że analizowane obiekty posiadają symetrię obrotową. Analizowany problem został przedstawiony w walcowym układzie współrzędnych (z, r, φ). Prądy w uzwojeniach płyną w kierunku wersora aφ . Rys. 2.2 przedstawia przekrój poprzeczny cewki wzdłuż płaszczyzny (r, z). Obszar Ω stanowi analizowaną cewkę, w której płynie prąd, zaś Γ jest obszarem z zadanym polem magnetycznym. W obszarze Ω prąd reprezentowany jest przez funkcję gęstości prądu (j = f (r, z)aφ , ∀(r, z) ⊆ Ω). Zadane prądy generują w obszarze Γ pole magnetyczne. Pole magnetyczne jest wektorem, który w rozważanym przypadku symetrii obrotowej posiada niezerowe składowe w kierunku r i z. W zagadnieniu kształtowania pola magnetycznego możemy wyodrębnić następujące podproblemy: X wyznaczenie rozkładu gęstości prądów j(z, r) w uzwojeniach cewki przy ustalonej powierzchni obrotowej Ω, X wyznaczenie kształtu powierzchni Ω przy ustalonym rozkładzie gęstości prądów w uzwojeniu,.

(13) 6. ROZDZIAŁ 2. PRZEDSTAWIENIE PROBLEMU. Rysunek 2.2: Przekrój poprzeczny wzdłuż płaszczyzny (r, z) obrazujący cewkę nawiniętą w obszarze Ω wytwarzającą pole magnetyczne w obszarze Γ.. X dobór zarówno kształtu jak i rozkładu gęstości prądu. Problemy powyższe można rozpatrywać nakładając ograniczenia na poszukiwane funkcje. Ograniczenia mogą dotyczyć zarówno kształtu uzwojenia cewki (np. Ω ⊆ [rcmin , rcmax ] × [zcmin , zcmax ]) jak i gęstości prądu (np. j ∈ [jcmin , jcmax ]).. 2.3. Pole magnetyczne cewki. W celu rozwiązania problemu syntezy konieczne jest (najczęściej wielokrotne) wyznaczenie pola magnetycznego pochodzącego od zadanego wymuszenia prądowego. Pole takie oblicza się w celu porównania go z polem zadanym. Do wyznaczenia pola magnetycznego można stosować metodę różnic skończonych, metodę elementów brzegowych [53] lub metodę elementów skończonych (MES)[9]. Podejście z wykorzystaniem MES jest najbardziej ogólne i pozwala wyznaczać pola elektromagnetyczne w dowolnych środowiskach. Podstawową wadą tej metody jest długi czas obliczeń. W pracy została zastosowana metoda podziału obszaru Ω na pojedyncze zwoje prądowe oraz zastosowanie wzorów analitycznych na wartość pola magnetycznego pochodzącego od przewodu kołowego. Podejście takie pozwala znacznie przyspieszyć obliczenia niezbędne do wyznaczenia pola w obszarze Γ. Metody tej nie można stosować w przypadku gdy pole magnetyczne wytworzone jest w środowisku nieliniowym. W dodatku A zostały wyprowadzone zależności określające wartości pola magnetycznego wytworzonego przez pojedynczy zwój kołowy (zależności (A.29), (A.30)). Przy założeniu, że analiza pola magnetycznego będzie wewnątrz cewki można założyć, że składowa pola magnetycznego w kierunku wersora z jest dużo.

(14) 7. 2.3. POLE MAGNETYCZNE CEWKI. większa od składowej pola magnetycznego w kierunku wersora r. Efekt ten jest znacznie silniejszy jeśli zadane pole magnetyczne jest polem jednorodnym. Zjawisko to określane jest jako “quadratic suppresion” [31]. Oczywiście w przypadku pola magnetycznego na osi z składowa Br wynosi zero. W niniejszej pracy rozważane będą problemy syntezy pola stacjonarnego dla obszaru Γ położonego na osi oraz poza nią. Analizowane będą również problemy syntezy pola niejednorodnego, ale tylko dla obszaru Γ położonego na osi z. Z dyskusji przedstawionej powyżej wynika, że w obu przypadkach wystarczy rozważać składową Bz wektora pola magnetycznego. W celu skrócenia zapisu składową Bz będziemy oznaczać symbolem B. Przy tych oznaczeniach na podstawie wzoru (A.29) pole magnetyczne B w punkcie o współrzędnych (r, z) wzbudzane przez prąd I pojedynczego zwoju o symetrii osiowej umiejscowionego w punkcie (r′ , z ′ ) jest równe:. ′. ′. B(r , z , r, z) =. √. (. µ0 I. 2π (r′ + rj )2 + (z ′ − z)2. ). (r′ )2 − r2 − (z ′ − z)2 K(k) + ′ E(k) (r − r)2 + (z ′ − z)2. (2.2) ′. r gdzie k 2 = (r′ +r)24r+(z ′ −z)2 zaś K(k) i E(k) są całkami eliptycznymi pierwszego i drugiego rodzaju. Przy założeniu, że pole indukowane jest na osi z układu współrzędnych (r = 0), równanie (2.2) redukuje się do postaci:. B(r′ , z ′ , 0, z) =. µ0 I(r′ )2 3. 2((r′ )2 + (z ′ − z)2 ) 2. .. (2.3). W celu analizy pola magnetycznego indukowanego w dowolnym punkcie obszaru Γ należy policzyć całkę po wszystkich pojedynczych zwojach znajdujących się w obszarze Ω. Wówczas kładąc I(r′ , z ′ ) = j(z ′ , r′ )dr′ dz ′ , gdzie j jest gęstością powierzchniową prądu w obszarze Ω, otrzymujemy: ∫∫. j(r′ , z ′ )K(r′ , z ′ ; r, z)dr′ dz ′. B(r, z) = Ω. (r′ , z ′ ) ⊆ Ω, (r, z) ⊆ Γ. (2.4). gdzie K(r′ , z ′ ; r, z) jest funkcją wiążącą geometrycznie rozważany zwój z prądem z punktem w którym wyznacza się pole magnetyczne. Równanie (2.4) pozwala obliczyć pole magnetyczne w dowolnym punkcie płaszczyzny (r, z). Rys. 2.3 przedstawia przekrój poprzeczny cewki. Cewka ta ograniczona jest przez płaszczyzny z = zmin oraz z = zmax oraz poprzez dwie krzywe obrotowe fl (z) oraz fu (z), które reprezentują ograniczenie od ”dołu” i od ”góry”. W rzeczywistych konstrukcjach najczęściej poszukuje się kształtu funkcji fl (z) ustalając górne ograniczenie jako stałe (fu (z) = const.). Wówczas korpus cewki otrzymuje zadany kształt natomiast cewka nawinięta jest równomiernie. W liniowej wersji problemu syntezy pola magnetycznego zadany jest obszar Ω oraz pole magnetyczne breq (r, z) w obszarze Γ i poszukiwany jest rozkład gęstości.

(15) 8. ROZDZIAŁ 2. PRZEDSTAWIENIE PROBLEMU. Rysunek 2.3: Cewka nawinięta na powierzchni obrotowej ograniczonej z dołu przez krzywą fl (z) i z góry przez fu (z). prądów w obszarze Ω. Wówczas równanie (2.4) jest niejednorodnym dwuwymiarowym liniowym równaniem całkowym Fredholma pierwszego rodzaju. Natomiast gdy poszukiwany jest kształt cewki wówczas równanie Fredholma pierwszego rodzaju traci swój liniowy charakter. Zagadnienie nieliniowe sprowadza się do znalezienia kształtu powierzchni obrotowej Ω. Problemy opisane równaniami Fredholma pierwszego rodzaju należą do klasy problemów źle postawionych. Głównymi metodami postępowania w przypadku takich problemów są metody regularyzacji, rozpowszechnione i rozpropagowane poprzez prace Tichonowa [35] [76] oraz Morozova [63]. Korzysta się również z innych metod, między innymi algorytmów stochastycznych [51]. Minimalizacja błędu średnio kwadratowego pomiędzy otrzymanym polem magnetycznym a polem zadanym jest równoważne minimalizacji funkcjonału: ∫∫ ( Γ. µ0 2π. )2. ∫∫ Ω. j(r′ , z ′ )K(r′ , z ′ , r, z)dz ′ dr′ − breq (r, z). dzdr. (2.5). Często w zagadnieniach syntezy pola magnetycznego rozważa się problem uzyskania stałego pola (breq = const.). Wówczas problem uzyskania zadanego pola magnetycznego w obszarze Γ jest równoważny wyznaczeniu jednorodnego pola magnetycznego na granicy ∂Γ obszaru Γ. Wynika to bezpośrednio z praw Gausa i Ampera (brak źródeł pola wewnątrz obszaru Γ). Równanie (2.5) opisuje przypadek ciągły problemu. Jedną z metod numerycznego rozwiązania problemu jest zastosowanie dyskretyzacji przestrzennej. W tym celu obszar Ω zostaje podzielony na n obszarów, zaś obszar Γ na m obszarów. Przykładowo, gdy obszar Ω jest prostokątem można go podzielić na n = w × k prostokątów. Dyskretyzację obszaru Γ uzyskuje się przez wybór m punktów w obszarze Γ lub w przypadku pola jednorodnego w zbiorze ∂Γ..

(16) 9. 2.4. PRZYPADEK LINIOWY. Dla uproszczenia obliczeń każdy z elementów podziału zbioru Ω można zastąpić pojedynczą cewką przez którą płynie prąd całkowity I. Zakładając, że obszar Ω podzielono na n pojedynczych elementów, natomiast obszar Γ na m elementów (najczęściej zakłada się m n) otrzymamy wersję dyskretną problemu opisanego równaniem (2.5). m ∑ j=1. (. n µ0 ∑ Ik K(rk , zk , rj , zj ) − breq (rj , zj ) 2π k=1. )2. (2.6). ∫∫. Można również całkę Ω j(r′ , z ′ )K(r′ , z ′ , r, z)dz ′ dr′ w równaniu (2.5) liczyć numerycznie zakładając podział obszaru Ω na pod obszary o stałej gęstości prądów. W szczególnych przypadkach przy wykorzystaniu symetrii można nawet wspomnianą całkę podwójną wyznaczyć analitycznie. Podejścia powyższe są równoważne. Minimalizując wyrażenie (2.6) otrzymamy rozwiązanie problemu dyskretnego. Podobnie jak w przypadku ciągłym możemy poszukiwać wartości prądów Ik , minimalizujących problem. Wówczas mamy do czynienia z problemem liniowym. Na poszukiwany rozkład prądów można nałożyć ograniczenia (np. I ∈ [imin , imax ]), wówczas problem staje się problemem liniowym z ograniczeniami. Możemy również poszukiwać kształtu obszaru Ω. Wtedy problem polega na znalezieniu wartości ri , zi i staje się problemem nieliniowym.. 2.4. Przypadek liniowy. W tym podrozdziale zostanie przedstawiony problem dla przypadku gdy poszukiwany jest rozkład prądów w cewce, natomiast położenia uzwojeń są ustalone. Załóżmy, że w obszarze Ω umieszczono n zwojów w punktach (rk , zk ), gdzie rk oznacza promień zwoju zaś zk jego współrzędną z. Załóżmy ponadto, że w obszarze Γ wybrano m > n punktów o współrzędnych (rj , zj ). Załóżmy, że wartość zadanego pola w punkcie o numerze j wynosi bj . Oznaczmy przez xk ∗ wartość natężenia prądu w zwoju o numerze k. Celem optymalizacji jest dobór wartości xk , tak aby indukowane pole magnetyczne było możliwie bliskie zadanemu polu bj w każdym z wybranych punktów obszaru Γ o współrzędnych (rj , zj ). Zgodnie z zależnościami (2.3) i (2.4) pole magnetyczne od pojedynczego zwoju zależy liniowo od prądu płynącego w tym zwoju. Można zatem sformułować następujące równanie liniowe: Ax = b (2.7) gdzie macierz A ∈ Rm×n jest macierzą współczynników, x ∈ Rn jest wektorem szukanych prądów, zaś wektor b ∈ Rm definuje zadane pole magnetyczne. Warto∗. ze względu na przejrzystość matematyczną standardowy symbol oznaczający prąd I został zamieniony przez symbol x.

(17) 10. ROZDZIAŁ 2. PRZEDSTAWIENIE PROBLEMU. ści współczynników macierzy A wynikają bezpośrednio z zależności (2.3) i (2.4). Ak,j =. (. µ0 xk. √. 2π (rk + rj )2 + (zk − zj )2. ). r2 − rj2 − (zk − zj )2 K(k) + k E(k) (rk − rj )2 + (zk − zj )2. (2.8) 4r r. k j gdzie k 2 = (rk +rj )2 +(z 2 . Jeśli nie nakładamy żadnych ograniczeń na wartości k −zj ) prądów xk to mamy do czynienia z problemem optymalizacji bez ograniczeń. W rzeczywistych problemach często w sposób naturalny pojawiają się jednak ograniczenia. Przykładowo zwykle żąda się, aby prądy płynęły w jednym kierunku (np. xk 0) lub żeby wartości prądów nie przekraczały pewnej wartości maksymalnej (|xk | ¬ xmax ). Prowadzi to do problemów z ograniczeniami.. 2.5. Błąd dyskretyzacji. W niniejszym podrozdziale dokonamy analizy błędu dyskretyzacji, tzn. różnicy pomiędzy polem pochodzącym od cewki punktowej oraz cewki prostokątnej ze stałą gęstością prądu. W przypadku pola magnetycznego generowanego w punkcie na osi z układu współrzędnych można wyznaczyć zależność analityczną. Wymiary cewki oraz obszar Γ są takie same jak użyte w późniejszych obliczeniach (tabela 4.1 w rozdziale 4.1).. 2.5.1. Cewka nieskończenie cienka. Rozważmy cewkę walcową nawiniętą nieskończenie cienkim przewodem. Promień cewki wynosi r. Zakładamy, że cewka jest nawinięta od punktu z1 do punktu z2 . W uzwojeniu tej cewki płynie prąd o stałej gęstości liniowej równej j = 1 A · m−1 . Wówczas na podstawie wzoru (2.3) pole magnetyczne w dowolnym punkcie z na osi współrzędnych możemy wyznaczyć z zależności: ∫. (. ). z − z1 z2 − z √ b(z) = +√ 2 2 2 r + (z − z ) r + (z − z2 )2 z1 1 (2.9) Dyskretyzacja przestrzeni cewki odpowiada zastąpieniu jej n zwojami. Prąd każdego zwoju wówczas wynosi: z2. µ0 j r2. µ0 j dξ = 3/2 2 2 2 2 (r + (z − ξ) ). x=j. z2 − z1 n. (2.10). W przypadku podziału na n zwojów wartość pola magnetycznego w dowolnym punkcie z na osi opisana jest zależnością: bn (z) =. n ∑. µ0 x r2. i=1. 2 (r2 + (z − ξi )2 )3/2. (2.11).

(18) 11. 2.5. BŁĄD DYSKRETYZACJI. 0. 10. r=0.20 m r=0.47 m r=0.73 m r=1.00 m. −5. δ. 10. −10. 10. −15. 10. 0. 1. 10. 2. 10. 10. 3. 10. n. Rysunek 2.4: Błąd dyskretyzacji dla cienkiej cewki walcowej dla różnych jej promieni. (. ). 1 gdzie ξi = z1 + z2 −z i − 12 dla i = 1, . . . , n jest pozycją zwoju o numerze i. n Poprzez b(z) będziemy rozumieli wektor określający wartość pola w dowolnym punkcie obszaru Γ. Załóżmy, że naszym obszarem Γ będzie odcinek usytuowany na osi z układu współrzędnych. Odcinek ten został podzielony na m = 1000 równomiernie rozłożonych punktów. Względny błąd dyskretyzacji definiujemy jako. δ=. ||b(z) − bn (z)||22 ||b(z)||22. (2.12). gdzie ∥ · ∥2 oznacza normę euklidesową. Na rys. 2.4 przedstawiono błąd dyskretyzacji w funkcji liczby pojedynczych zwojów. Wartości pola wyznaczone były dla m = 1000 punktów na osi z. Wykres został narysowany w skali logarytmicznej. Łatwo zauważyć, że wraz ze wzrostem liczby zwojów błąd ten szybko maleje i dla n = 10 błąd jest mniejszy niż 10−5 . W związku z powyższym wartość pola reprezentowana przez pojedyncze uzwojenia jest obarczona dość małym błędem już dla niewielkiej liczby uzwojeń. Z wykresu można również wywnioskować, że im większy jest promień cewki tym błąd jest mniejszy.. 2.5.2. Cewka o skończonej grubości. Podobnie można przeprowadzić analizę w przypadku, gdy cewka nie jest idealnie cienka, a posiada pewną grubość. W tym wypadku przekrój poprzeczny cewki stanowi prostokąt Ω = [z1 , z2 ] × [r1 , r2 ]. Zakłada się, że w cewce płynie prąd o stałej gęstości powierzchniowej j = 1 A · m−2 . Wówczas pole na osi wyraża się.

(19) 12. ROZDZIAŁ 2. PRZEDSTAWIENIE PROBLEMU. 0. 10. −2. 10. −4. 10. δ. −6. 10. −8. 10. nr=1 nr=126. −10. 10. nr=251 nr=375 n =500. −12. r. 10. 0. 1. 10. 2. 10. 10. 3. 10. n. z. Rysunek 2.5: Błąd dyskretyzacji cewki walcowej w funkcji podziału w kierunku z przy ustalonym podziale w kierunku r.. wzorem ∫. r2. ∫. z2. b(z) = r1. z1. µ0 jr2. µ0 j dξdr = 3/2 2 2 (r2 + (z − ξ)2 ). (. Ψ1 z log − z1 log Ψ1 + z2 log Ψ2 Ψ2. ). (2.13) √ r2 +(z−zi )2 +r2 dla i = 1, 2. Na rys. 2.5 i 2.6 przedstawiono błąd gdzie Ψi = √ 22 2 r1 +(z−zi ) +r1. dyskretyzacji dla podziału w kierunkach z i r† . Z przeprowadzonej analizy wynika, że już dla podziału cewki na stosunkowo niewielką liczbą elementów błąd dyskretyzacji jest niewielki i w obliczeniach praktycznych może zostać pominięty.. 2.6. Definicja funkcji celu. W calu sformułowania problemu syntezy pola jako zadania optymalizacyjnego należy zdefiniować funkcje celu. Najczęściej funkcję celu definiuje się w sposób następujący: f (x) = ∥b(x) − breq ∥ †. iz. (2.14). gdzie nr i nz oznaczają podział cewki na pojedyncze zwoje odpowiednio w kierunku osi r.

(20) 13. 2.6. DEFINICJA FUNKCJI CELU. 0. 10. −2. 10. −4. 10. n =1 z. nz=126. δ. −6. 10. nz=251 nz=375 nz=500. −8. 10. −10. 10. −12. 10. 0. 1. 10. 2. 10. 10. 3. 10. n. r. Rysunek 2.6: Błąd dyskretyzacji cewki walcowej w funkcji podziału w kierunku r przy ustalonym podziale w kierunku z.. gdzie ∥ · ∥ jest pewną normą zadaną w przestrzeni Rm . Zwykle stosuje się jedną z następujących norm: ||y||1 =. m ∑. |yi | (norma suma modułów). (2.15). i=1. v um u∑ ||y||2 = t yi2. (norma euklidesowa). (2.16). (norma maksimum). (2.17). i=1. ||y||∞ = max |yi | i. Rozwiązanie problemu optymalizacji polega na znalezieniu minimum globalnego funkcji celu. W przypadku stosowania metod gradientowych wygodne jest użycie normy euklidesowej ∥ · ∥2 z uwagi na jej gładkość. W przypadku normy euklidesowej często w celu uniknięcia nieróżniczkowalnej operacji pierwiastkowania jako funkcję celu stosuje się równoważne z punktu widzenia problemu optymalizacji wyrażenie, f (x) = ∥b(x) − breq ∥22. (2.18). Tak zdefiniowana funkcja celu będzie poddana optymalizacji w całej niniejszej pracy. W przypadku optymalizacji liniowej b(x) = Ax, gdzie macierz A jest opisana wzorem (2.8). Zatem w przypadku liniowym problem optymalizacji staje się zagadnieniem najmniejszych kwadratów..

(21) 14. 2.7. ROZDZIAŁ 2. PRZEDSTAWIENIE PROBLEMU. Ocenianie wyników. W celu porównania wyników obliczeń niezbędne jest wprowadzenie wielkości dodatkowo oceniających rezultat obliczeń. Wyniki obliczeń można oceniać na dwa sposoby. Pierwszy polega na porównaniu rozkładu różnicy pomiędzy polem magnetycznym otrzymanym a polem zadanym. Ocenę taka może spełniać wektor błędu względnego ε = [ε1 , ε2 , . . . , εm ] zdefiniowany wzorem: εi = |(bi − bi,req )/bi,req | ,. i ∈ [1, . . . , m]. (2.19). gdzie bi jest otrzymaną wartością pola w i – tym punkcie obszaru Γ, natomiast bi,req oznacza zadaną wartość pola w tym punkcie. Wektor ε zawiera pełną informację o błędzie. Ze względów praktycznych wygodniej jest porównywać wielkości skalarne. W związku z tym w niniejszej pracy będziemy dokonywać porównania wyników na podstawie następujących współczynników: γ1 = ∥ϵ∥2 ∥ϵ∥1 γ2 = m γ3 = ∥ϵ∥∞. (2.20) (2.21) (2.22). W przypadku liniowym i gdy pole zadane jest polem jednorodnym współ√ czynnik γ1 = f (ˆ x), gdzie x ˆ jest wyznaczonym rozkładem prądów w cewce. Współczynnik γ1 określa wartość normy euklidesowej błędu rozwiązania. Natomiast współczynnik γ2 określa wartość średnią błędu względnego pomiędzy polem otrzymanym a polem zadanym. Z kolei współczynnik γ3 określa maksymalny błąd względny rozwiązania. W calu porównania własności rozwiązań problemu liniowego będziemy używać następujących wielkości: β1 = ∥ˆ x∥22 =. n ∑. x ˆ2i .. (2.23). β2 = ∥ˆ x∥∞ = max |ˆ xi |.. (2.24). i=1 i∈1,...,n. Natomiast w przypadku problemu optymalizacji kształtu cewki będziemy używać współczynnika: ∫ j 2 dΩ. β1 =. (2.25). Ω. Ze względu na fakt, że w przypadku optymalizacji kształtu cewki wartość gęstości prądu w całym obszarze Ω pozostaje niezmienna współczynnik β2 jest stały i porównywanie go nie ma sensu. Współczynnik β1 jest proporcjonalny do wartości energii czy to zgromadzonej w uzwojeniu cewki czy to traconej w jej uzwojeniach. W związku z tym, będziemy go nazywali współczynnikiem energetycznym. Współczynnik β2 jest równy maksymalnej wartości prądu rozwiązania..

(22) 2.7. OCENIANIE WYNIKÓW. 15. Posiada on również ważne uzasadnienie konstrukcyjne związane z ewentualnymi naprężeniami mechanicznymi występującymi w uzwojeniu, szczególnie istotnymi w przypadku konstrukcji cewek opartych o uzwojenia wykonane z nadprzewodnika. Należy podkreślić, że przedstawione współczynniki mają charakter pomocniczy i są wykorzystywane wyłącznie do oceny jakości rozwiązań. Nie podlegają one procesowi optymalizacji. W niniejszej pracy nie stosuje się optymalizacji wielokryterialnej, która brałaby więcej elementów pod uwagę..

(23) Rozdział 3. Metody optymalizacji. W rozdziale niniejszym przedstawione zostaną wybrane metody optymalizacji, które mogą być zastosowane do rozwiązania problemu syntezy pola magnetycznego.. 3.1. Podstawowe definicje. Celem optymalizacji jest znalezienie wektora x ∈ Rn , dla którego funkcja celu f (x) : D → R osiąga ekstremum w pewnym ustalonym zbiorze D ⊆ Rn . Jeśli D = Rn to mamy do czynienia z optymalizacją bez ograniczeń, w przeciwnym wypadku mówimy o optymalizacji z ograniczeniami. Zbiór D może być zarówno ograniczony jak i nieograniczony. Każdy punkt ze zbioru D nazywamy punktem dopuszczalnym. Zbiór D często definiowany jest za pomocą pewnego zbioru ograniczeń. Mogą to być zarówno ograniczenia równościowe (gi (x) = 0) jak i nierównościowe (gi (x) 0). W dowolnym punkcie x ∈ D ograniczenia te mogą być podzielone na dwa rozłączne podzbiory. Pierwszy zawiera ograniczenia równościowe – nazywane aktywnymi, zaś drugi ograniczenia nierównościowe – zwane wolnymi. Jeśli jakieś aktywne ograniczenie jest zbyteczne, tzn. otrzymujemy takie same rozwiązanie bez niego, to rozwiązanie takie nazywamy zdegenerowanym. Rozpoczniemy od przypomnienia definicji wybranych pojęć z dziedziny optymalizacji używanych w dalszej części pracy. Niech D będzie dowolnym podzbiorem w przestrzeni Rn , zaś f dowolną funkcją o wartościach rzeczywistych zdefiniowana na zbiorze D. f : Rw ⊃ D → R Będziemy rozważać problem minimalizacji funkcji f na zbiorze D, gdzie: I. D = Rn – problem optymalizacji bez ograniczeń. 16.

(24) 3.2. OPTYMALIZACJA Z OGRANICZENIAMI. 17. II. D = {x ∈ Rn : g1 (x) = 0, . . . , gm (x) = 0}, gdzie g1 , . . . , gm : Rn → R – problem optymalizacji z ograniczeniami równościowymi. III. D = {x ∈ Rn : g1 (x) ¬ 0, . . . , gm (x) ¬ 0}, gdzie g1 , . . . , gm : Rn → R – problem optymalizacji z ograniczeniami nierównościowymi. W naszym przypadku interesujące problemy zawarte są w podpunktach I oraz III. W szczególności najczęściej będziemy rozważać ograniczenia kostkowe typu xmin ¬ xk ¬ xmax . Punkt x0 ∈ D nazywamy minimum globalnym funkcji f na zbiorze D jeśli: ∀x ∈ D. f (x) f (x0 ).. Punkt x0 ∈ D nazywamy minimum lokalnym funkcji f jeśli: ∃ε > 0 ∀x ∈ D. ∥x − x0 ∥ < ϵ ⇒ f (x) f (x0 ).. Minimum nazywamy ścisłym, jeżeli w powyższych definicjach zachodzi f (x) > f (x0 ), dla x ̸= x0 . Jeśli zamienimy funkcję f funkcją −f wówczas problem minimalizacji zamieni się w problem maksymalizacji. W związku z tym analizuje się zwykle tylko problemy minimalizacji. Niepusty podzbiór D ⊂ Rn nazywamy wypukłym jeśli dla każdego x1 , x2 ∈ D zachodzi {λx1 + (1 − λ)x2 : λ ∈ [0, 1]} ⊂ D. Funkcję f nazywamy wypukłą na zbiorze D jeśli dla każdego x1 , x2 ∈ D oraz λ ∈ [0, 1] spełniona jest nierówność f (λx1 + (1 − λ)x2 ) ¬ λf (x1 ) + (1 − λ)f (x2 ). Funkcja f jest ściśle wypukła jeśli nierówność jest ostra dla λ ∈ (0, 1). Jeżeli mamy do czynienia z ograniczeniami liniowymi, to zbiór ograniczeń jest zbiorem wypukłym. Jeżeli mamy do czynienia z minimalizacją funkcji ściśle wypukłej z wypukłym zbiorem ograniczeń to mamy do czynienia z tzw. programowaniem wypukłym. W takim przypadku istnieje minimum i jest on jedyne. Punkt x ˆ spełniający zależność ∇f (ˆ x) = 0 nazywamy punktem krytycznym. [. ]T. (x) , . . . , ∂f Wielkość ∇f (x) = ∂f∂x(x) jest gradientem funkcji f w punkcie x. ∂xn 1 W przypadku optymalizacji bez ograniczeń, przy założeniu, że funkcja celu jest ciągła i dwukrotnie różniczkowalna w punkcie x ˆ, warunkiem koniecznym aby punkt x ˆ był lokalnym ścisłym minimum jest spełnienie warunku ∇f (ˆ x) = 0 oraz dodatnia określoność macierzy H(ˆ x). Macierz H(ˆ x) jest Hessianem funkcji 2f f , tzn. macierzą drugich pochodnych cząstkowych, Hi,j = ∂x∂i ∂x . Dodatkowo j jeśli funkcja f jest ściśle wypukła to punkt x ˆ spełniający powyższe warunki jest minimum globalnym funkcji f .. 3.2. Optymalizacja z ograniczeniami. Rozwiązywanie zagadnień optymalizacji związanych z rzeczywistymi problemami często wymaga uwzględniania różnorodnych ograniczeń. W niniejszej pracy.

(25) 18. ROZDZIAŁ 3. METODY OPTYMALIZACJI. ograniczeniami są objęte zarówno wartości prądów w uzwojeniach oraz położenie uzwojeń związane z ograniczeniami konstrukcyjnymi. Uwzględnianie ograniczeń polega na takiej modyfikacji funkcji celu, aby jej minimalizacja prowadziła do rozwiązania problemu z ograniczeniami. Zdefiniujmy problem optymalizacji z ograniczeniami, gdzie minimalizuje się funkcję f (x), x ∈ Rn z p ograniczeniami równościowymi hj (x) oraz z m ograniczeniami nierównościowymi gi (x). x ˆ = arg min. (3.1). f (x). x∈Rn. gi (ˆ x) ¬ 0,. i = 1, . . . , m. hj (ˆ x) = 0,. j = 1, . . . , p. Częstą metodą uwzględniania ograniczeń w funkcji celu jest metoda Lagrange’a. W metodzie tej wprowadza się tak zwaną funkcję Lagrange’a, która dla problemu (3.1) ma postać: L(x, λ, µ) = f (x) +. m ∑ i=1. λi gi (x) +. p ∑. µj hj (x). (3.2). j=1. Liczby λi , µj nazywamy współczynnikami Lagrange’a związanymi z ograniczeniami nierównościowymi i równościowymi. Twierdzenie Lagrange’a mówi, że jeśli funkcja L ma ciągłe pochodne cząstkowe, a punkt x ˆ jest rozwiązaniem optyˆ µ malnym problemu z ograniczeniami, to istnieją wektory λ, ˆ takie, że wszystkie ˆ pochodne cząstkowe funkcji Lagrange’a L w punkcie (ˆ x, λ, µ ˆ) są równe zeru. Rozwinięciem twierdzenie Lagrange’a jest twierdzenie Karusha-Kuhna-Tuckera [52], które formułuje warunki istnienia globalnego minimum dla problemu programowania wypukłego, oraz pozwala na rozwiązanie zagadnienia optymalizacji z ograniczeniami. Twierdzenie to mówi, że punkt x ˆ jest globalnym minimum problemu programowania wypukłego jeśli istnieje λ = (λ1 , . . . , λm )T ∈ Rm oraz µ = (µ1 , . . . , µp )T ∈ Rp takie, że: 1. ∇x L(ˆ x, λ, µ) = ∇f (ˆ x) + λT ∇g(ˆ x) + µT ∇h(ˆ x) = 0 2. gi (ˆ x) ¬ 0,. i = 1, . . . , m. 3. hj (ˆ x) = 0,. j = 1, . . . , p. 4. λi 0,. i = 1, . . . , m. 5. λi gi (ˆ x) = 0,. i = 1, . . . , m. Powyższe warunki noszą nazwę warunków Karusha-Kuhna-Tuckera (KKT). Inną metodą uwzględniania ograniczeń w problemach optymalizacyjnych jest metoda funkcji kary. Rozróżnia się metodę kary wewnętrznej i zewnętrznej. Metody polegają na skonstruowaniu dodatkowej funkcji kary S : Rn → R+ , której wartości zależą od analizowanego punktu. Jeśli punkt jest punktem dopuszczalnym.

(26) 3.3. KLASYFIKACJA METOD OPTYMALIZACJI. 19. (mieści się w zbiorze ograniczeń) wówczas funkcja ta przyjmuje małe wartości, natomiast kiedy analizowany punkt nie spełnia warunków ograniczeń wówczas jej wartość znacznie wzrasta. Dodając wspomnianą funkcję do funkcji celu problemu otrzymamy nowy problem optymalizacyjny bez ograniczeń. F (X) = f (x) + ξS(x) gdzie F (X) jest nową funkcją celu, f (x) jest funkcją celu problemu bez ograniczeń, zaś S(x) jest funkcją kary uwzględnioną ze współczynnikiem wagi ξ. W przypadku zewnętrznej funkcji kary wartość jej wynosi zero, gdy punkt x jest punktem dopuszczalnym, natomiast jest dodatnia w przypadku przeciwnym. Dla wewnętrznej funkcji kary jej wartość jest większa lub równa zero dla punktów dopuszczalnych, natomiast jej wartość znacznie wzrasta gdy punkt zbliża się do granicy ograniczeń. Zarówno metoda Lagrange’a jak i metoda funkcji kary powodują zastąpienie problemu z ograniczeniami problemem bez ograniczeń. Metoda Lagrange’a wraz z warunkami KKT jest dużo bardziej odpowiednia w przypadku programowania wypukłego i sprawdza się w przypadku zastosowania metod gradientowych optymalizacji. Natomiast metody funkcji kary są często stosowane w przypadku algorytmów niedeterministycznych.. 3.3. Klasyfikacja metod optymalizacji. W rozdziale tym zostaną przedstawione metody optymalizacji zastosowane w niniejszej pracy. Klasyfikacja metod optymalizacji jest ściśle związana z rodzajem problemu, w szczególności zależy od tego, czy problem jest liniowy czy nieliniowy, czy problem jest z ograniczeniami czy bez ograniczeń. Rozpatrzmy problem minimalizacji bez ograniczeń w postaci: znaleźć takie x ˆ aby f (ˆ x) = minn f (x) x∈R. W metodach tych można wyróżnić następujące przypadki: 1. Dysponujemy jedynie algorytmem wyznaczającym wartość funkcji celu f (x), oraz spodziewamy się, że funkcja f (x) ma bardzo nieregularny charakter i w związku z tym obliczenie gradientu ∇f (x) obarczone będzie dużym błędem. Wówczas należy skorzystać metod bezgradientowych. 2. Dysponujemy zarówno algorytmem wyliczającym wartość funkcji celu f (x) oraz algorytmem wyznaczania jej gradientu ∇f (x). Oznacza to, że dysponujemy formułą analityczną wyliczającą pochodne cząstkowe ∂f∂x(x) dla i i = 1, . . . , n. Wówczas zalecane są metody gradientowe..

(27) 20. ROZDZIAŁ 3. METODY OPTYMALIZACJI. 3. Dysponujemy jedynie algorytmem wyznaczania wartości funkcji f (x), ale wartość gradientu ∇f (x), choć nie jest on dany w postaci formuł analitycznych, może być z dużą dokładnością przybliżana przez ilorazy różnicowe. Postępowanie takie nosi nazwę estymacji gradientu. Możliwość estymacji gradientu występuje wtedy, gdy wartości funkcji f są obliczane z dużą dokładnością. Celowe jest wówczas rozwiązanie tego zadania za pomocą metod gradientowych, wykorzystujących metodę estymacji gradientu. Ważny jest również podział metod ze względu na kształt funkcji celu. Możemy mieć tutaj do czynienia z problemami liniowymi, które prowadzą przy minimalizacji normy euklidesowej residuum do zagadnienia najmniejszych kwadratów, oraz z problemami nieliniowymi. Wśród problemów nieliniowych ważnym podzbiorem jest podzbiór problemów gdy funkcja celu jest funkcją wypukłą. Podsumowując, metody optymalizacji możemy podzielić na metody gradientowe, w których wykorzystuje się nie tylko wartości funkcji celu w danym punkcie poszukiwań ale również wartości jej gradientu czy nawet hessjanu. Natomiast w metodach bezgradientowych korzysta się tylko i wyłącznie z wyliczonej wartości funkcji celu. Ważną podgrupą metod bezgradientowych są metody losowe, które w trakcie poszukiwań minimum funkcji wykorzystują czynnik losowy. Zadanie optymalizacji staje się bardziej skomplikowane w przypadku analizy funkcji, która ma więcej niż jedno minimum lokalne. Wówczas stosowanie metod gradientowych może okazać się nieefektywne, ponieważ algorytmy te korzystają z lokalnych własności funkcji. Często otrzymujemy w rozwiązaniu ekstremum lokalne, które nie musi być poszukiwanym ekstremum globalnym. Wówczas można skorzystać z metody wielostartowej, która jest połączeniem metody gradientowej i losowej. Metody losowe charakteryzują się tym, że wyznaczenie rozwiązania zagadnienia optymalizacyjnego odbywa się w wyniku iteracji, polegających na przeszukiwaniu otoczenia aktualnego (w k – tej iteracji) punktu przybliżenia x(k) . Sposób tego przeszukiwania zależy od rodzaju metody. Często metody te są mało efektywne, ale mogą być stosowane wówczas, gdy metody gradientowe zawodzą. Jako głównych reprezentantów metod gradientowych można wymienić następujące metody [36]: X największego spadku i jej modyfikacje, X gradientów sprzężonych, X Newtona X quasi-Newtona Natomiast pośród wielu bezgradientowych wymienia się następujące metody: X Hooka i Jeevesa,.

(28) 3.4. ZAGADNIENIE NAJMNIEJSZYCH KWADRATÓW. 21. X Rosenbrocka, X sympleksową Neldera i Meada, X Powella i jej modyfikację, Do najprostszych metod losowych można zaliczyć metodę Monte Carlo, w której dokonuje się wielu losowań punktu z przestrzeni potencjalnych rozwiązań. Jako rozwiązanie optymalne przyjmuje się najlepsze znalezione rozwiązanie spośród wszystkich wylosowanych. Inną metodą losową jest algorytm błądzenia przypadkowego. Bieżące rozwiązanie jest oceniane i na podstawie jego położenia losuje się następny punkt, w taki sposób, że punkt bieżący jest wartością oczekiwaną rozkładu prawdopodobieństwa wykorzystywanego do generowania kolejnego punktu. Przez cały proces poszukiwań musi być pamiętane najlepsze z dotychczas znalezionych rozwiązań. Natomiast metody niedeterministyczne wyróżniają się zwykle dwiema cechami: pierwszą z nich jest prostota, drugą zaś duża efektywność działania w przestrzeniach i problemach, w których klasyczne metody optymalizacji i poszukiwania, nie dają zadowalających wyników. Często problemem może być nieciągłość funkcji celu, oraniczenie rozwiązania do zbioru liczb całkowitych itp. Spośród wielu algorytmów niedeterministycznych można wymienić następujące: X symulowane wyżarzanie, X algorytmy genetyczne i ewolucyjne, X algorytmy mrówkowe, X algorytmy oparte o zastosowanie sztucznych sieci neuronowych.. 3.4. Zagadnienie najmniejszych kwadratów. W nienijszym podrozdziale zostanie omówione zagadnienie najmniejszych kwadratów, gdyż pojawia się ono naturalnie przy sformułowaniu problemu syntezy pola magnetycznego przy użyciu normy euklidesowej w definicji funkcji celu. Zaprezentowane zostanie zagadnienie najmniejszych kwadratów bez ograniczeń, jak i z nałożonymi ograniczeniami.. 3.4.1. Zagadnienie najmniejszych kwadratów bez ograniczeń. W podrozdziale 2.4 została przedstawiona liniowa zależność pomiędzy poszukiwanymi prądami a zadanym polem magnetycznym (2.7). Ponieważ macierz współczynników A jest macierzą prostokątną, to mamy do czynienia z równaniem nadokreślonym. Wówczas o rozwiązaniu problemu mówi się w sensie najmniejszych kwadratów. Zagadnienie najmniejszych kwadratów polega na znalezieniu.

(29) 22. ROZDZIAŁ 3. METODY OPTYMALIZACJI. wektora x ˆ, który minimalizuje sumę kwadratów elementów wektora residualnego r(x) = b − Ax: x ˆ = arg min ||b − Ax||2. (3.3). x∈Rn. Dla problemu najmniejszych kwadratów funkcja celu ma postać f (x) = ||b − Ax||22 =. m ∑. (. bj −. j=1. n ∑. )2. Aj,i xi. .. (3.4). i=1. Funkcja ta osiąga minimum w punkcie w którym wartość jej pochodnej wynosi zero. . . n ∑ ∂f (x) = 2 bj − Aj,i xi  (−Ai,j ) = 0 ∂xi j=1. dla. i = 1, 2, . . . , n. (3.5). Po przekształceniach otrzymujemy tzw. układ równań normalnych m ∑ n ∑. Aj,i Aj,k xk =. j=1 k=1. m ∑. Aj,i bj ,. dla. i = 1, 2, . . . , n. (3.6). j=1. W zapisie macierzowym układ równań normalnych ma postać (. ). AT A x = AT b. (3.7). Równanie (3.7) jest równaniem z kwadratową macierzą AT A. Równanie to można rozwiązać bezpośrednio x = A† b, gdzie A† = (AT A)−1 AT nazywa się macierzą pseudoodwrotną (macierzą Moore’a-Penrose’a). Równanie normalne dla problemu najmniejszych kwadratów jest równaniem liniowym z macierzą kwadratową. Metody rozwiązywania takich równań dzieli się na: X Metody bezpośrednie – eliminacja Gaussa (i jej różne warianty) – metody dekompozycji (Choleskiego, QR, SVD) X Metody iteracyjne – metody iteracji prostej: Jackobiego, Gaussa-Seidla, wersje relaksacyjne – metody dla macierzy rzadkich. Trudność rozwiązania układu równań liniowych zależna jest od tzw. uwarunkowania problemu, którą ocenia się obliczając współczynnik uwarunkowania problemu. Współczynnik uwarunkowania dla macierzy kwadratowej A definiuje się jako κ(A) := ∥A−1 ∥2 ∥A∥2 =. maxi |λi | mini |λi |. (3.8).

(30) 3.4. ZAGADNIENIE NAJMNIEJSZYCH KWADRATÓW. 23. gdzie λi są wartościami własnymi macierzy A. Natomiast w przypadku, gdy macierz A jest macierzą prostokątną współczynnik uwarunkowania κ definiuje się jako maxi |σi | κ(A) := ||A† ||2 ||A||2 = (3.9) mini |σi | gdzie σ oznacza wartości szczególne macierzy prostokątnej A. Wartości szczególne definiuje się jako pierwiastki wartości własnych kwadratowej macierzy AT A. Jeśli wartość współczynnika uwarunkowania jest wysoka mówimy wówczas o złym uwarunkowaniu macierzy. Wówczas małe zaburzenia danych wejściowych powodują duże różnice w rozwiązaniu. W przypadku rozważanych problemów nie ma konieczności stosowania metod iteracyjnych, które najczęściej znajdują zastosowanie dla dużych układów równań lub macierzy rzadkich. Bezpośrednie rozwiązanie układu równań normalnych Jedną z metod rozwiązania zagadnienia najmniejszych kwadratów jest bezpośrednie rozwiązanie układu równań normalnych t.j. obliczenie AT A, a następnie rozwiązania równania liniowego z tą macierzą. Można to zrobić za pomocą metody eliminacji Gaussa. Metoda bezpośrednia dobrze działa w przypadku gdy wskaźnik uwarunkowania macierzy AT A jest mały. Ponieważ AT A jest macierzą symetryczną można również zastosować rozkład Choleskiego. Wówczas procedura rozwiązania problemu najmniejszych kwadratów ma postać: C = AT A, d = AT b, RT R = C, RT y = d, Rx = y. gdzie trójkątną macierz R oblicza się za pomocą rozkładu Choleskiego. Wynik obliczeń zależy od stosowanej reprezentacji zmiennoprzecinkowej liczb. W standardzie IEEE 754 liczby rzeczywiste podwójnej precyzji przedstawione są z błędem względnym nie większym niż u = 21 β t−1 , β = 2, t = 53 [42]. Taką reprezentację stosowano w obliczeniach opisanych w niniejszej pracy. Rozwiązanie problemu najmniejszych kwadratów metodą dekompozycji Choleskiego staje się niestabilne w przypadku gdy wartość wskaźnika uwarunkowania κ(A) jest bliska u−0.5 (dla standardu IEEE754 w wersji podwójnej precyzji ∼ 108 )[2]. W przypadku gdy wskaźnik uwarunkowania macierzy κ(AT A) jest duży, znacznie lepsze efekty można uzyskać stosując rozkład QR macierzy A do rozwiązania układu równań normalnych. Dekompozycja macierzy QR Metoda dekompozycji QR macierzy polega na rozkładzie macierzy A na iloczyn dwóch macierzy Q i R. A = QR, (3.10) gdzie macierz Q ∈ Rm×n jest macierzą ortogonalną (QT Q = In ), natomiast macierz R ∈ Rn×n jest macierzą trójkątną górną. W celu dokonania dekompozycji.

(31) 24. ROZDZIAŁ 3. METODY OPTYMALIZACJI. macierzy QR najczęściej stosuje się metody wykorzystujące ortogonalizację Grama–Schmidta, obroty Gibsa lub odbicia Hausholdera [47][77]. W naszym przypadku metoda Hausholdera jest bardziej stabilna. Dysponując rozkładem A = QR rozwiązanie układu równań normalnych można uzyskać rozwiązując układ z macierzą trójkątną Rx = QT b. Można to zauważyć dokonując przekształcenia układu równań normalnych. AT Ax = AT b, (QR)T QRx = (QR)T b, RT QT QRx = RT QT b, Rx = QT b. Istotnym elementem tej metody jest fakt uniknięcia wyliczania źle uwarunkowanej macierzy AT A, co znacznie poprawia stabilność numeryczną problemu w porównaniu do metody opisanej w poprzednim podrozdziale. Mimo większego nakładu pracy potrzebnej do przeprowadzenia dekompozycji QR macierzy A, wyniki są bardziej stabilne i w efekcie metoda może być stosowana dla większych wymiarów. Dekompozycja SVD Do analizy problemu zastosowano również rozkład macierzy A względem wartości szczególnych (SVD – singular value decomposition) [47]. Rozkład ten polega na przedstawieniu macierzy A ∈ Rm×n w postaci iloczynu: A = U ΣV T. (3.11). gdzie Σ = diag(σ1 , . . . , σr , 0, . . . , 0) ∈ Rm×n , σ1 . . . σr 1, r jest rzędem macierzy A, natomiast macierze U ∈ Rm×m , V ∈ Rn×n są ortonormalne (U T U = V T V = I oraz ∥U ∥2 = ∥V ∥2 = 1). Po dokonaniu dekompozycji można rozwiązać równanie normalne rozwiązując równanie: ΣV T x = U T b. 3.4.2. (3.12). Zagadnienie najmniejszych kwadratów z ograniczeniami. Nieujemne zadanie najmniejszych kwadratów Celem nieujemnego zadania najmniejszych kwadratów jest wyznaczenie wektora x ∈ Rn o nieujemnych elementach minimalizującego funkcję celu f (x) = ||Ax − b||22 dla zadanych A ∈ Rm×n , b ∈ Rm , co możemy zapisać jako: x ˆ = arg min ∥Ax − b∥22 ,. (3.13). x∈Rn , x0. W literaturze angielskojęzycznej powyższy problem określany jest jako “nonnegative least square problem” (NNLS). Warunki KKT dla powyższego problemu.

(32) 3.4. ZAGADNIENIE NAJMNIEJSZYCH KWADRATÓW. 25. mają postać x 0, ∇f (x) = AT (Ax − b) 0, (. )T. ∇f (x)T x = AT (Ax − b). (3.14). x = 0.. W związku z powyższym problem sprowadza się do wyznaczenia nieujemnego wektora x spełniającego równanie (Ax − b)T Ax = 0. Metody rozwiązujące nieujemne zadanie najmniejszych kwadratów można podzielić na algorytmy z aktywnym zbiorem ograniczeń (zastosowane w niniejszej pracy), metody gradientowe (projected quasi-Newton NNLSQ, metoda Landwebera) oraz inne (metoda punktu wewnętrznego, principal block pivoting method). W przypadku ograniczeń nieujemnych zbiór ograniczeń ma postać g(x) = {g1 (x), . . . , gm (x)}, gdzie gi (x) = −x. Mówimy, że i−te ograniczenie jest aktywne w danym punkcie dopuszczalnym x jeśli gi (x) = 0. W przeciwnym wypadku ograniczenie to nazywamy pasywnym. Zbiór indeksów ograniczeń aktywnych w danym punkcie x oznaczamy przez A(x). Pozostałe indeksy związane z ograniczeniami pasywnymi pozostają w zbiorze P(x). Zbiór P jest dopełnieniem A, tzn. A(x) ∪ P(x) = {1, . . . , m}. Algorytmy oparte na aktywnym zbiorze ograniczeń polegają na rozwiązaniu problemu najmniejszych kwadratów ograniczonego tylko do zmiennych, dla których ograniczenia są pasywne. Pozostałe zmienne przyjmują wartość ograniczeń równościowych. Algorytm rozwiązujący nieujemne zadanie najmniejszych kwadratów metodą aktywnego zbioru ograniczeń został zaproponowany przez Lawsona i Hansona R [54]. Jest on zaimplementowany w środowisku Matlab ⃝ jako funkcja lsqnonneg [73]. W algorytmie tym rozwiązanie konstruowane jest iteracyjnie, przy czym algorytm jest zbieżny w skończonej liczbie kroków. Algorytm 1 przedstawia klasyczny algorytm NNLS (Non-negative Least Squares). Macierz AP oznacza zredukowaną macierz A ograniczoną tylko do kolumn odpowiadających zbiorowi pasywnemu P. Lawson i Henson udowodnili, że algorytm znajduje jednoznaczne rozwiązanie po skończonej liczbie iteracji. Niestety nie istnieją metody pozwalające oszacować liczbę iteracji potrzebną do rozwiązania problemu. Czasem warto algorytm przerwać ponieważ im algorytm jest bliżej rozwiązania tym poprawa jest wolniejsza. Podczas działania algorytmu układ równań normalnych [(AP )T AP ]sP = (AP )T b (linia 10 algorytmu) rozwiązuje się metodą bezpośrednią. Może to powodować problemy ze stabilnością numeryczną, szczególnie dla macierzy źle uwarunkowanych. W celu uniknięcia tych problemów do rozwiązania układu równań normalnych zastosowano metodę rozkładu QR. Po tej zmianie linia 10 algorytmu (1) ma postać: AP = QP RP – dekompozycja QR macierzy AP sP = (RP )−1 (QP )T b.

(33) 26. ROZDZIAŁ 3. METODY OPTYMALIZACJI. Algorytm 1 Algorytm NNLS (Lawson i Hanson). Input: A ∈ Rm×n , b ∈ Rm , ε Output: x ˆ = arg minx0 ∥Ax − b∥2 1: P = ∅, A = {1, 2, . . . , n}, x = 0, w = AT (b − Ax) 2: while A ̸= ∅ and maxi∈A (wi ) > ε do 3: j = arg maxi∈A (wi ) 4: P = P ∪ {j}, A = A \ {j} 5: sP = [(AP )T AP ]−1 (AP )T b 6: while min(sP ) ¬ 0 do i 7: α = − mini∈P xix−s i 8: x := x + α(s − x) 9: Wybierz elementy A oraz P 10: sP = [(AP )T AP ]−1 (AP )T b 11: sR = 0 12: end while 13: x=s 14: w = AT (b − Ax) 15: end while. Dekompozycja pozwala uniknąć niestabilnych numerycznie operacji kosztem wydłużenia czasu obliczeń. Rozkład macierzy AP stosowany był tylko w wypadku jej dużego wymiaru. Algorytm 2 Algorytm FNNLS (Bro i de Jung). Input: A ∈ Rm×n , b ∈ Rm Output: x ˆ = arg minx0 ∥Ax − b∥2 1: P = ∅, A = {1, 2, . . . , n}, x = 0, w = AT b − (AT A)x 2: while A ̸= ∅ and maxi∈A (wi ) > tolerance do 3: j = arg maxi∈A (wi ) 4: P = P ∪ {j}, A = A \ {j} 5: sP = [(AT A)P ]−1 (AT b)P ) 6: while min(sP ) ¬ 0 do i 7: α = − mini∈P xix−s i 8: x := x + α(s − x) 9: Wybierz elementy A oraz P 10: sP = [(AT A)P ]−1 (AT b)P ) 11: sA = 0 12: end while 13: x=s 14: w = AT (b − Ax) 15: end while Bro i Jong [10] zaprezentowali udoskonalony algorytm nazwany Fast NNLS (FNNLS). Algorytm 2 jest modyfikacją klasycznego algorytmu NNLS zapropono-.

(34) 27. 3.4. ZAGADNIENIE NAJMNIEJSZYCH KWADRATÓW. waną przez Bro i de Junga [10]. Idea polega na wykorzystaniu w każdym kroku algorytmu raz obliczonych pełnych macierzy AT A, AT b. Wersja ta powoduje skrócenie czasu obliczeń szczególnie dla problemów o dużym wymiarze. Z uwagi na brak możliwości połączenia tej wersji z metodą rozkładu QR w niniejszej pracy nie była ona stosowana. Następnie pojawiło się jeszcze ulepszenie zwany fast combinatorial NNLS [7] również skracające czas działania algorytmu poprzez pewne modyfikacje w obliczeniach. Dokładniejszy opis tych metod można znaleźć w [8] [10] [13] [48] [49]. Metoda najmniejszych kwadratów z ograniczeniami kostkowymi W poprzednim podrozdziale na rozwiązanie problemu najmniejszych kwadratów nałożono ograniczenie, aby było ono nieujemne. Inne ograniczenie, często pojawiające się w zastosowaniach praktycznych, polega na tym, że wynik ma należeć do ustalonego przedziału. Takie ograniczenie nazywać będziemy ograniczeniem kostkowym. Wówczas problem najmniejszych kwadratów z ograniczeniami ma postać: x ˆ = arg min ∥Ax − b∥22 (3.15) x∈Rn , xi ∈[li ,ui ]. gdzie −∞ < li < ui < ∞, i = 1, . . . , n, A ∈ Rm×n , b ∈ Rm . Ponieważ ||Ax − b||22 = (Ax − b)T (Ax − b) = xT AT Ax − 2bT Ax + bT b (. 1 T T 1 =2 x A Ax − (AT b)T x + bT b 2 2. ). to wprowadzając oznaczenia Q = AT A ∈ Rn×n oraz q = −AT b ∈ Rn , problem (3.15) można przekształcić do postaci: x ˆ=. 1 T x Qx + q T x x∈Rn ,xi ∈[li ,ui ] 2 arg min. (3.16). Równanie powyższe opisuje problemem programowania kwadratowego z ograniczeniami kostkowymi [3] [18] [32] [61] [78] [79]. Macierz Q jest macierzą kwadratową dodatnio półokreśloną. Jeśli, co zwykle ma miejsce, macierz Q jest macierzą dodatnio określoną, to przy liniowych ograniczeniach mamy zagwarantowane istnienie i jednoznaczność rozwiązania (na podstawie własności programowania wypukłego opisanej w rozdziale 3.1). Ograniczenia kostkowe można przedstawić jako ograniczenia nierównościowe opisane za pomocą liniowej funkcji G(x) = Cx − d gdzie macierz C ∈ R2n×2n jest macierzą blokową diagonalną cii = 1 dla i ∈ 1, . . . , n, cii = −1 dla i ∈ w + 1, . . . , 2n natomiast wektor d określa górne i dolne ograniczenia: [ ] di = [ li , di+w ] = −ui dla i ∈ {1, . . . , n}. Korzystając z I 0 l oznaczeń C = ,d= problem (3.16) można opisać jako problem 0 −I −u.

(35) 28. ROZDZIAŁ 3. METODY OPTYMALIZACJI. programowania kwadratowego z ograniczeniami nierównościowymi: 1 x ˆ = arg min xT Qx + q T x, x∈Rn 2 Cx d.. (3.17). Do rozwiązania powyższych zagadnień, w przypadku niedużych wymiarów problemu często stosuje się metody oparte na aktywnym zbiorze (active set methods). Natomiast stosuje się również inne metody takie jak: metody punktu wewnętrznego (interior-point methods), projekcyjne metody gradientowe (projected gradient methods). Podobnie jak przy problemie nieujemnym zagadnieniu najmniejszych kwadratów, w niniejszej pracy użyto metod opartych na aktywnym zbiorze. Dla problemu (3.17) funkcja Lagrange’a ma postać: 1 L(x, λ, µ) = xT Qx + xT q − λT (x − l) + µT (x − u) (3.18) 2 Gdzie λ, µ są współczynnikami Lagrange’a reprezentującymi ograniczenia dolˆ µ ne i górne. Warunki KKT (Karusha-Kuhna-Tuckera) konieczne aby x ˆ, λ, ˆ były rozwiązaniem problemu, mają postać: ˆ+µ Qˆ x+q−λ ˆ=0 ˆ λ 0, µ ˆ0 ˆ x − l) = 0 λ(ˆ. (3.19). µ ˆ(u − x ˆ) = 0 l¬x ˆ¬u Wykorzystując warunki (3.19) można zastosować metodę aktywnego zbioru do rozwiązania problemu. Algorytm 3 przedstawia procedurę rozwiązania problemu najmniejszych kwadratów przy ograniczeniach kostkowych [78]. W odróżnieniu od algorytmów rozwiązujących nieujemne zadanie najmniejszych kwadratów mamy do czynienia z dwoma zbiorami aktywnymi (L i U). Zbiory te zawierają indeksy niewiadomych x, które nie mieszczą się w zbiorze ograniczeń (odpowiednio dolnym i górnym). Wówczas nadaje im się wartości graniczne (odpowiednio li lub ui ) i poszukuje rozwiązania operując tylko na zbiorze pasywnym S. W każdej iteracji zachodzi związek L∪U ∪S = {1, . . . , n}. Algorytm kończy działanie, gdy wektory x, µ, λ spełniają warunki KKT, czyli są rozwiązaniem problemu (linia 13 algorytmu). Zakończenie działania następuje również po przekroczeniu maksymalnej liczby iteracji kmax . W linii 9 przedstawionego algorytmu, rozwiązywany jest układ równań liniowych: n ∑ j=1. (k+1). Qi,j xk+1 + qi = λi j. (k+1). − µi. ,. dla i ∈ {1, . . . , n}. (3.20).

(36) 29. 3.4. ZAGADNIENIE NAJMNIEJSZYCH KWADRATÓW. Algorytm 3 Algorytm aktywnego zbioru dla problemu programowania kwadratowego z ograniczeniami kostkowymi Input: Q ∈ Rn×n , q ∈ Rn , l ∈ Rn , u ∈ Rn , kmax Output: x ˆ = minx 21 xT Qx + q T x oraz li ¬ xi ¬ ui , ∀i ∈ {1, . . . , n} 1: k = 0, λ(0) = 0, µ(0) = 0, x(0) = −Q−1 q. 2: for k = 0{to kmax do } (k) (k) (k) 3: L(k) = i : xi < li lub {xi = li oraz λi 0} {. (k). (k). U (k) = i : xi. 5:. S (k) = {1, 2, . . . , n}\(L(k) ∪ U (k) ). 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18:. {. (k+1). > ui lub {xi. (k). 4:. = ui oraz µi }. 0}. }. (k+1). xi = li , µi = 0, ∀i ∈ L(k) (k+1) (k+1) xi = ui , λi = 0, ∀i ∈ U (k) (k+1) (k+1) λi = 0, µi = 0, ∀i ∈ S (k) (k+1) (k+1) (k+1) rozwiąż Qx + q = λi − µi ze względu na n niewiadomych (k+1) (k) xi , ∀i ∈ S (k+1) λi , ∀i ∈ L(k) (k+1) µi , ∀i ∈ U (k) (k+1) (k+1) if li ¬ xi ¬ ui , ∀i ∈ S (k) oraz µi 0, ∀i ∈ U (k) oraz λi 0, ∀i ∈ (k) L then STOP x ˆ = x(k+1) . else k ←k+1 end if end for (k+1). (k+1). Ponieważ dla i ∈ S (k) mamy λi = µi = 0, równanie (3.20) możemy łatwo (k+1) (k) rozwiązać ze względu na xi , dla i ∈ S biorąc pod uwagę aktywne zbiory określone w liniach 3 i 4 algorytmu: ∑. (k+1). Qi,j xj. =−. j∈S (k). ∑. ∑. Qi,j lj −. j∈L(k). Qi,j uj − qi ,. dla i ∈ S (k). (3.21). (k+1). , dla i ∈. j∈U (k) (k+1). Natomiast współczynniki Lagrange’a λi , dla i ∈ L(k) oraz µi U (k) można wyznaczyć bezpośrednio z następujących wzorów: (k+1). λi. =. ∑. (k+1). + qi , dla i ∈ L(k). (3.22). (k+1). − qi , dla i ∈ U (k). (3.23). Qi,j xj. j∈{1,...,n} (k+1). µi. =−. ∑. Qi,j xj. j∈{1,...,n}. Opisaną powyżej metodę zastosowano do rozwiązania problemów z ograniczeniami kostkowymi rozważanych w niniejszej pracy..

(37) 30. ROZDZIAŁ 3. METODY OPTYMALIZACJI. 3.5. Metoda regularyzacji. Popularną metodą rozwiązywania źle uwarunkowanych równań liniowych jest metoda regularyzacji [24][72]. Metoda ta pozwala otrzymać bardziej stabilne i przede wszystkim gładkie rozwiązania. Często przy zwiększaniu wymiaru rozpatrywanego problemu obserwowane są następujące własności macierzy A: 1. wartości własne macierzy A zanikają stopniowo do zera, oraz pojawia się coraz więcej wartości własnych bliskich zeru, 2. stosunek największej wartości własnej do najmniejszej wartości własnej jest bardzo duży. Świadczą one o złym uwarunkowaniu problemu, czyli o dużej czułości rozwiązania na wartość wektora b. Małe zmiany wartości wektora b (np. błąd reprezentacji zmiennoprzecinkowej) powodują duże różnice w rozwiązaniu. Metoda regularyzacji jest metodą rozwiązania problemu polegającą na zmianie funkcji celu, tak aby uzyskać wymagane własności rozwiązania. W odróżnieniu od metod podanych poprzednio głównym celem jest otrzymanie bardziej gładkich i stabilnych rozwiązań. Najbardziej znaną metodą regularyzacji jest regularyzacja Tichonowa [76][63]. Ideą tej metody jest modyfikacja funkcji celu przez dodanie członu regularyzującego. W przypadku zadania najmniejszych kwadratów zmodyfikowana funkcja celu ma postać: (. ). xλ = arg min ∥Ax − b∥22 + λ2 ∥L(x − x∗ )∥22 ,. (3.24). x. gdzie L jest macierzą Tichonowa (najczęściej jest to macierz jednostkowa) kontrolującą gładkość rozwiązania, x∗ jest przypuszczalnym (oczekiwanym) rozwiązaniem problemu, natomiast parametr regularyzacyjny λ jest współczynnikiem kontrolującym wpływ regularyzacji na rozwiązanie. Jeśli współczynnik λ jest duży to informacja o oryginalnym problemie staje się mniej istotna, natomiast małe λ oznacza, że składnik regularyzacyjny ma mniejsze znaczenie. Przy założeniu L = I oraz x∗ = 0, stosując podobne przekształcenia jak w przypadku metody najmniejszych kwadratów, otrzymujemy układ równań normalnych: ( ) AT A + λ2 I x = AT b (3.25) Powyższe równanie najczęściej rozwiązuje się za pomocą rozkładu SVD macierzy A lub za pomocą metod iteracyjnych [40]. Głównym problemem jest wyznaczenie optymalnej wartości parametru regularyzacji λ. Spośród wielu metod wyznaczania parametru λ, jedną z najczęściej stosowanych metod jest wykorzystanie wykresu przedstawiającego normę ∥Lxλ ∥2 w funkcji normy ∥Axλ − b∥2 . Rys. 3.1 przedstawia przykładowy wykres takiej krzywej. Krzywa ta często nazywana jest krzywą L (ang. ”L-curve”) ze względu na jej kształt przypominający.

(38) 31. 3.5. METODA REGULARYZACJI. log x. 2. wiksze , rozwizanie bardziej gadkie. mniejsze , mniejszy wpyw regularyzacji. log A x-b. 2. Rysunek 3.1: Przykładowy wykres dla metody krzywej L wyboru wartości parametru λ.. literę L. Optymalnego parametru λ oczekuje się w narożniku krzywej. Miejsce to jest kompromisem pomiędzy oczekiwanymi własnościami rozwiązania, a błędem rozwiązania problemu wprowadzonym przez obecność składnika regularyzującego. Jak zostanie pokazane w dalszej części pracy zastosowanie metody krzywej L w odniesieniu do rozważanych problemów nie przynosi zadowalających wyników. Jest to przede wszystkim związane z faktem, że konstruowane wykresy nie mają typowego kształtu krzywej L. W związku z tym zaproponowano inną metodę doboru optymalnej wartości parametru regularyzacyjnego λ. Jako parametr optymalny uznano taką najmniejszą wartość λ, dla której rozwiązanie problemu (3.24) jest nieujemne dla każdej współrzędnej, co można formalnie opisać jako: (. ). λopt = min{λ : arg min ∥Ax − b∥22 + λ2 ∥x∥22 0}.. (3.26). x. Zastosowanie tej metody powoduje automatyczne spełnienie warunku, że otrzymane rozwiązanie jest nieujemne (tak jak w nieujemnej metodzie najmniejszych kwadratów). Obecność członu regularyzacyjnego zwiększa gładkość i stabilność rozwiązania. W celu wyznaczenia wartości optymalnej λopt zastosowano metodę bisekcji (Algorytm 4). Algorytm ten dla ustalonej macierzy A oraz wektora b znajduje optymalną wartość współczynnika regularyzacji λopt z założonym z góry błędem ε. Algorytm rozpoczyna poszukiwanie optymalnego współczynnika od zadanej wartości początkowej λ0 (zakłada się, że λ0 > λopt , tzn. rozwiązanie uzyskane dla λ0 jest nieujemne)..