[WikiDyd] [TitleIndex] [WordIndex

Laboratorium Metod Numerycznych II
Zakład Elektrotechniki Teoretycznej i Informatyki Stosowanej
Instytut Elektrotechniki Teoretycznej i Systemów Informacyjno Pomiarowych, PW

Ćwiczenie Nr 2

ROZWIĄZANIE ZAGADNIENIA ODWROTNEGO PRZY WYKORZYSTANIU MES

Autor: R. Szmurło

Cel ćwiczenia

Rozwiązanie zadania odwrotnego na przykładzie lokalizacji obiektu znajdującego się w środowisku przewodzącym na podstawie pomiarów potencjałów zdjętych z kilkunastu elektrod. Ćwiczenie jest przykładem tomografii impedancyjnej. Oto zdjęcie rzeczywistego urządzenia z laboratorium podstaw elektromagnetyzmu.

attachment:wanna.png.jpg

Uwaga! Do wykonania ćwiczenia niezbędny jest działający solver MES, który został wykonany w ćwiczeniu 1.

Literatura

  1. T. Guziak, A. Kamińska, B. Pańczyk, J. Sikora, “Metody Numeryczne w Elektrotechnice”, Rozdział 10: Rozwiązywanie równań różniczkowych cząstkowych
  2. S. Bolkowski, i inni, “Komputerowe metody analizy pola elektromagnetycznego”, Rozdział 4: Metoda elementów skończonych (MES)

Wstęp teoretyczny

Rozwiązanie prezentowanego zagadnienia odwrotnego polega na minimalizacji pewnej funkcji celu. W naszym przypadku funkcja celu powinna osiągnąc minimum wtedy gdy badany obiekt będzie znajdował się w oryginalnej pozycji.

W ćwiczeniu zakładamy, że znamy wartości funkcji w 12 ściśle określonych węzłach siatki. (Czyli otrzymujemy je od prowadzącego zajęcia.) Węzły te reprezentują elektrody pomiarowe. Zadaniem jest odnalezienie takiego położenia obiektu wewnątrz, dla którego w wyniku przeprowadzonej przez nas symulacji, uzyskamy możliwie najbliższe wartości funkcji w wspomnianych 12 węzłach. Dla tak postawionego problemu możemy sformułować funkcję celu następująco:

/usr/local/bin/latex['--interaction=nonstopmode', 'latex_5c9bd1fcf6d82966f25c331099e5739e181d184d_p.tex']/tmp/tmpZDz7Jtlatex error! exitcode was 3 (signal 0), transscript follows:

[Wed Jul 17 06:24:59 2013] [error] failed to exec() /usr/local/bin/latex

gdzie x i y to położenie badanego obiektu,

/usr/local/bin/latex['--interaction=nonstopmode', 'latex_b24ae4dc0d3e2f18bb651adecc1d303a17240b1f_p.tex']/tmp/tmpfxR0U1latex error! exitcode was 3 (signal 0), transscript follows:

[Wed Jul 17 06:25:01 2013] [error] failed to exec() /usr/local/bin/latex
to wartości funkcji w węzłach podane prze prowadzącego (czyli zadane) oraz
/usr/local/bin/latex['--interaction=nonstopmode', 'latex_ee23d72aed21b63b59e405278e7316e701c73d6a_p.tex']/tmp/tmpdWTU5Olatex error! exitcode was 3 (signal 0), transscript follows:

[Wed Jul 17 06:25:03 2013] [error] failed to exec() /usr/local/bin/latex
to wartości funkcji, które uzyskaliśmy po wykonaniu naszej symulacji dla konkretnej pozycji obiektu x,y.

Oto schematyczny obrazek (na razie 24 elektrody, ale zobaczymy w wersji ostatecznej ile ich będzie):

attachment:schemat2.png

Zakres ćwiczenia

Do wykonania ćwiczenia potrzebne będzie:

Zadanie 1.

Rozwiązać zadanie dostarczone przez prowadzącego (dostarczone macierze N i E) za pomocą własnego symulatora. (1 punkt)

Kod programu zwracającego macierze węzłów i elementów: attachment:prostokat2.m

Wywołanie:

function [N,E] = prostokat2(xPos,yPos)
% xPos moze byc w zakresie 0.25 do 5.75
% yPos moze byc w zakresie 0.25 do 3.75
%
% N - bedzie zawierac trzy kolumny:
%     1 - wspolrzedne x
%     2 - wspolrzedne y
%     3 - znacznik warunku brzegowego:
%        0 - Neumann = 0 (czyli nic nie robimy :-)
%        1 - Dirichlet D1 (lewa elektroda)
%        2 - Dirichlet D2 (prawa elektroda)
%
% E - bedzie zawierac 4 kolumny:
%     1-3 - indeksy wezlow danego elementu
%     4 - wartosc wspolczynnika materialowego c z rownania Laplace'a

dstep = 0.25; % okresla gestosc generowanej siatki
c1 = 1;       % wspolczynnik materialowy zewnetrznego obszaru
c2 = 100;     % wspolczynnik mater. poszukiwanego obszaru

Podpowiedź: Skorzystaj z poniższego schematu kodu, uzupełniając odpowiednie fragmenty.

function zadanie2(px,py)

[N,E] = prostokat2(px,py);

(... tutaj jest rozwiazanie MES ...)

hold on
scatter3(N(1:14,1),N(1:14,2),x(1:14), 'filled');
trisurf(E(:,1:3), N(:,1), N(:,2), x,x);
hold off

% okresl szukane wartosci funkcji w wybranych 14 pierwszych wezlach
xc = [
   (....)
    ];
x(1:14)

% wartosc fukcji celu:
dot(x(1:14)-xc,x(1:14)-xc)

Zadanie 2.

Wyznaczyć wrażliwość funkcji celu w trzech wskazanych przez prowadzącego punktach. (1 punkt)

Podpowiedź: Wrażliwość funkcji celu wyznaczmy obliczając gradient funkcji celu w kierunku osi x oraz osi y. Gradient w kierunku osi x wyznaczamy wykonując dwie symulacje: dla określonej pozycji x=x1 oraz y=y1, oraz po przesunięciu obiektu o pewien DeltaX, czyli; x=x1+DeltaX oraz y=y1. Oznacza to, że otrzymuje dwie wartości funkcji celu: f(x1,y1) = F1 oraz f(x1+DeltaX, y1) = F2. Wrażliwość funkcji jest wyznaczana z wzoru: W = (F2-F1) / DeltaX

Zadanie 3.

Napisać program, który znajdzie lokalizację obiektu na podstawie wartości funkcji ui podanych przez prowadzącego wykorzystując metodę Newtona. (1 punkt)

Przykładowe wartości funkcji ui (xc - to x celu):

xc = [
    0.8321
    1.3354
    1.8104
    3.9836
    6.6366
    8.3052
    9.1367
    9.1588
    8.3518
    6.7577
    4.3807
    2.3065
    1.1517
    0.5875
    ];

Uwaga! W trakcie zajęć będzie trzeba wykonać jedno z poniższych zadań. Które zadanie będzie obowiązywać na konkretnych zajęciach będzie podane na początku danych zajęć.

Zadanie 4.

Napisać program wykorzystujący algorytm genetyczny podany prze prowadzącego do znalezienia rzeczywistej lokalizacji. (2 punkty)

Podpowiedzi:

Wykonując zadanie należy zaimplementować wskazany przez prowadzącego algorytm selekcji, krzyżowania lub mutacji. Należy zwrócić szczególną uwagę na pliki ga_rs.m, który zawiera w komentarzach informacje o zastosowanym kodowaniu, plik funk_celu.m, który definiuje funkcję używaną w procedurze oceny populacji oraz funkcje reprodukcja.m oraz krzyzowanie_i_mutacja.m. Duża część informacji zostanie podana w trakcie zajęć.

Uwaga! Aby zadanie nie liczyło się zbyt długo, należy zmniejszyć dyskretyzację zadania MES. Należy wyedytować plik prostokat2.m i odpowiednia zwiększyć dyskretyzację np.  dstep = 0.5  .


2015-09-23 06:44