Zastosowanie identyfikacji opartej na wnioskowaniu rozmytym do identyfikacji modeli dynamicznych na przykładzie modelu protezy POLVAD

W pracy został przedstawiony sposób modelowania dynamicznych właściwości procesu z wykorzystaniem parametrycznych modeli rozmytych. Podejście to zostało zaprezentowane na przykładzie badania dynamiki sztucznej protezy serca POLVAD, przeprowadzonego na podstawie licznych eksperymentów testowych. We wstępie opisano badany proces i podkreślono pewne właściwości o charakterze nieciągłego działania, których występowanie ogólnie ogranicza zastosowanie opisu transmitancyjnego. Wyprowadzono podstawowe zależności do prowadzenia obliczeń parametrycznych modeli dynamicznych drogą oszacowań statystycznych dla modeli liniowych parametrycznych oraz modeli rozmytych. Przedstawiono wyniki badań optymalnego modelu liniowego i modelu rozmytego oraz przeprowadzono obszerną dyskusję uzyskanych modeli. Uzyskany opis dynamiki badanego procesu odtwarza wszystkie istotne zachowania badanego procesu w sposób całkowicie wystarczający do prowadzenia badań analizujących różne algorytmy sterowania lub regulacji, jak również umożliwia wprowadzenie np. optymalizacji pracy zaworów sterujących.

Modelowanie dynamiki procesów staje się jednym z podstawowych sposobów badania układów w fazie tworzenia i badania prototypowych układów oraz ich sterowania. Technika ta uzyskała już swoją nazwę – fast prototyping – i jest podstawowym narzędziem dla projektowania nowoczesnych układów oraz badania rozwiązań alternatywnych, poprzedzających wykonanie docelowych konstrukcji. Podstawą tego podejścia jest utworzenie zbioru modeli odtwarzających zachowanie komponentów układu. Modele takie można opracowywać na drodze analitycznych rozważań lub posługiwać się techniką identyfikacji. Jakkolwiek opis uzyskany drogą identyfikacji statystycznej wykazuje pewne ograniczenia, np. nie można z jego pomocą dokonywać optymalizacji działania zespołu, którego działanie zastępuje, ale ma również określone zalety, np. w przypadku istniejącego, gotowego zespołu może bardzo wiernie modelować jego zachowanie w układzie [2]. Przykład takiego zastosowania identyfikacji (opracowanie modelu dynamiki sztucznej protezy serca POLVAD) zostanie przedstawiony w artykule.

Modelowany proces – proteza POLVAD

Proteza POLVAD jest mechanicznym zespołem, który czasowo zastępuje pracę wybranej komory serca. Napęd tego urządzenia jest pneumatyczny, a zasada działania oparta na odtworzeniu cyklicznego, pulsacyjnego charakteru pompowania krwi do układu krążenia. Jej konstrukcję opracowano w Polsce i jest obecnie wykorzystywana w wielu szpitalach. Fotografia samej komory została pokazana na rys. 1.

Na rys. 2 są przedstawione oznaczenia podstawowych wielkości, umożliwiających właściwy opis tego zespołu oraz wielkości fizyczne mierzone podczas eksperymentów. Zespół POLVAD zasilany jest przepływem F_a powietrza o niskim poziomie ciśnienia P_zas, które mechanicznie przepycha elastyczną przeponę S_k i powoduje wyrzut krwi F_wyl przez otwartą zastawkę wylotową z_wyl, gdy ciśnienie krwi P_k jest wyższe od ciśnienia krwi w kaniuli wylotowej P_wyl. W przypadku cyklu zasysania krwi ciśnienie P_wl powinno być większe od P_k, co spowoduje otwarcie zastawki wlotowej z_wl i przepływ F_wl krwi do komory POLVAD. Oczywiście odpowiednie warunki zapewnia właściwe wysterowanie ciśnieniowe od strony zawory pneumatycznego P_a, który prowadzi do ciągłej wymiany strumienia powietrza F_a. Określony sposób sterowania czasowego powinien zapewniać całkowite opróżnienie komory z krwi (V_k=0 zapewnia brak zalegania krwinek i zapobiega zakrzepom) oraz odpowiednie opróżnienie komory powietrznej V_a w zależności od zadanego poziomu wspomagania protezy POLVAD, a określonej przez zadaną objętość martwą komory powietrznej V . Opisany zespół zawiera bezwładnościowo działające zastawki mechaniczne, które zamykają lub otwierają określone kanały w sposób zbliżony do działania diod i tak były do tej pory modelowane [1, 5]. Zasada konstrukcji oraz działanie poszczególnych elementów wydają się proste i łatwe do zamodelowania mechanicznego, jednak obserwacje i przeprowadzone pomiary (P_wyl, F_wyl, P_zas, F_wl, P_wl) ujawniają zjawiska, które nie są zgodne z tak prostym podejściem do zachodzących zjawisk.

Na rys 3. przedstawiono przykładowo pomierzone (co 1 ms) zmiany ciśnień P_k, P_zas oraz przepływu F_wyl. Początek fazy wyrzutu jest zgodny z oczekiwaniami: P_zas przewyższa P_wyl (obie wielkości wyskalowane odrębnie!) i następuje wyrzut krwi F_wyl > 0 (dla lepszej rozróżnialności poziom 0 został zaznaczony dodatkową przerywaną linią).

Następnie poziom ciśnienia zasilającego maleje, przepływ wylotowy F_wyl zmniejsza się do dużej ujemnej wartości, a następnie wzrasta (zaznaczono te fazy przerywanymi strzałkami). Oznacza to, że zastawki nie działają jak diody oraz, że należy uwzględnić dodatkowe zjawiska przepływowe związane ze spójnością strugi, dynamiką zmian ciśnień, nieznaną zmienną bezwładnością strugi krwi itd.

W takim przypadku modelowanie zjawisk dynamicznych przestaje być prostym zadaniem, tym bardziej że badany zespół POLVAD mógł pracować w dość szerokim przedziale zmian: zadawanej średniej różnicy ciśnień, prędkości zmian (czasu cyklu) i procentowego czasu wypełnienia. W prowadzonych badaniach dobry model przepływu wylotowego był jednym z podstawowych czynników warunkujących prawidłowe zaprojektowanie całego układu sterowania, który powinien zapewniać odpowiednie warunki zasilania obwodu krążenia, przy jednoczesnej minimalizacji zużycia energii w celu zwiększenia czasu resursu urządzenia podczas zasilania autonomicznego. Poszukiwany model powinien uzależniać wielkość wyrzutu F_wyl od warunków zewnętrznych: ciśnienia wlotowego P_wl i przepływu wlotowego F_wl, ciśnienia wylotowego P_wyl oraz oczywiście sygnału sterującego P_zas (rys. 4).

Wielkości ΔP oraz τ oznaczają odpowiednio zadaną wartość różnicy ciśnień (średnie wylotowe względem średniego wlotowego) oraz czas wypełnienia cyklu – długość fazy wyrzutowej. Pozyskanie dobrego opisu, w formie podanej na rys. 4, zdecydowanie poprawiłoby jakość prowadzonych prac projektowych. Obok technik modelowych [1] zdecydowano się wykorzystać technikę identyfikacji rozmytej.

Rozmyty model dynamiki procesu

Załóżmy, że proces o przedstawionej wyżej strukturze można lokalnie przybliżyć za pomocą zestawu transmitancji operatorowych G_Fwyl Vj(s), zatem dla zerowych warunków początkowych sygnał wyjściowy z powyższego procesu byłby opisany zależnością:

(1)

gdzie transmitancje operatorowe G_Fwyl Vj(s) są określone przez odpowiednie współczynniki oraz wartości opóźnień:

(2)

i nieznane są wartości współczynników typu α oraz β jak również wartości opóźnień T_j. Wyznaczenie wartości tych współczynników wykonywane jest najczęściej w sposób pośredni metodą identyfikacji parametrycznej, tzn. w początkowym etapie wyznaczane są transmitancje dyskretne [2]:

(3)

na podstawie danych pomierzonych z określonym krokiem próbkowania Δ w taki sposób, aby modelowany sygnał wyjściowy był jak najbliższy pomierzonym wartościom F_wyl(kΔ), które uzależnione są od pomierzonych wartości sygnałów wejściowych P_wyl(kΔ), P_wl(kΔ), F_wl(kΔ), P_zas(kΔ) oraz nieznanych parametrów b_0j, ..., a_j:

(4)

gdzie tzw. wektor wejść do modelu v(kΔ) należący do zbioru R^L zawiera sygnały pomierzone w odpowiednich (z uwzględnieniem opóźnień) chwilach czasowych (dla uproszczenia dalszych zapisów symbol Δ został pominięty)

(5)

a wektor kolumnowy θ należący do zbioru R^L zawiera wartości nieznanych współczynników transmitancji (3)

(6)

Wektor współczynników θ może być oszacowany różnymi metodami. Jedną z najprostszych jest metoda najmniejszej sumy (LS) kwadratów błędów modelu (4), która pozwala dokonać oszacowania na podstawie znajomości kolumnowego wektora pomierzonych wyjść Y = [F_wyl(1),…, F_wyl(N)]^T należący do zbioru R^N, gdzie N oznacza liczbę pomiarów wielkości wyjściowej, oraz macierzy V należący do zbioru R^NxL zawierającej N wierszowych wektorów wejść do modelu (4) w odpowiednich chwilach czasowych V^T=[v(1)^T, …, v(N)^T]

(7)

Zaletą tej metody jest prosty, nierekursywny algorytm obliczeń. Następnie wartości współczynników wektora θ_LS są przekształcane drogą aproksymacji na współczynniki transmitancji G_{Fwyl Vj}(s) (2).

Model wyznaczony w taki sposób może bardzo dokładnie określać właściwości dynamiczne badanego procesu [3, 4], ale ma właściwości lokalne i w większości przypadków wymaga wyznaczenia wielu modeli, np. dla różnych warunków zasilania ΔP lub różnych czasów wypełnienia. Jednocześnie, w przypadku gdy występują zjawiska nieciągłe, np. działanie zaworu odcinającego, wówczas modele liniowe o opisie (2) lub (3) stają się niedokładne.

Bardzo korzystne w takim przypadku jest zastosowanie modeli rozmytych [1, 2, 6], tzn. modeli o takiej samej strukturze jak w (4), ale dla których wektor współczynników θ jest zmienny i uzależniony od tzw. zmiennej r rozmywającej (lub kilku zmiennych r₁, r₂, ...). Zmienna rozmywająca jest określona przez funkcje przynależności φ_i, (rys. 5). Na przykład, w przypadku, gdy zadana dla protezy POLVAD średnia różnica ciśnień ΔP między ciśnieniem na wylocie P_wyl oraz na wlocie P_wl będzie z obszaru pomiędzy 0 a 100 mmHg, wówczas obowiązywać będzie model θ₁ (wektor współczynników) stosowny dla małej różnicy ciśnień - wyznaczony dla partycji θ₁. Podobnie będzie w przypadku, gdy zmienna rozmywająca ΔP przyjmie wartość 145 mmHg, wówczas należy zastosować model θ₂ wyznaczony dla średniej różnicy ciśnień ΔP, odpowiadający funkcji przynależności φ₂.

Natomiast w przypadku zmiennej rozmywającej ΔP równej ~110 mmHg funkcja przynależności φ₁ jest równa np. 0,7, a funkcja φ₂jest równa 0,3 i wówczas model wynikowy jest kombinacją dwóch modeli θ = φ₁θ₁ + φ₂θ₂ = 0,7·θ₁ + 0,3·θ₂.

Ogólnie można wyznaczyć model rozmyty [2] w formie kombinacji liniowej modeli wyznaczonych dla różnych przedziałów funkcji przynależności. Przykładowo, dla dwóch zmiennych (w funkcji dyskretnej chwili czasu k) rozmywających r₁(k) i r₂(k), o odpowiednio p i q przedziałach funkcji przynależności postać modeli będzie następująca

(8)

Funkcje przynależności muszą spełniać określone warunki przy zastosowaniu powyższego formalizmu

(9)

dla dowolnej wartości zmiennej rozmywającej r(k). Model taki nazywany jest modelem TSK (Takagi-Sugeno-Kang).

Algorytm oszacowania wektorów modeli cząstkowych qi, przy wykorzystaniu modeli rozmytych, jest bardziej rozbudowany niż formuła (7) i wymaga jednoczesnego obliczenia wektorów dla wszystkich modeli θ_j, tzn. wyznaczenia wektora θ^T= [θ^T₁, θ^T₂, …, θ^T_L] złożonego z L wektorów θ_i, który zapewni minimum wartości wskaźnika

(10)

Powyższe zagadnienie ma rozwiązanie w innej postaci [6]

(11)

i wymaga odwracania znacznie większych macierzy niż oszacowanie (7). Z tego powodu, przy tej formie modelu ważne jest tzw. uwarunkowanie numeryczne obliczeń, które określa poprawne wyznaczanie modeli cząstkowych. Wyznaczenie rodziny modeli cząstkowych (11) nie jest złożone obliczeniowo, a podstawową trudność stanowi właściwy dobór kształtu liczby funkcji przynależności. Istnieją również algorytmy [4, 7], które umożliwiają automatyczny dobór funkcji przynależności φ_i w taki sposób, by różnica między pomierzonym sygnałem y i jego oszacowaniem na podstawie modelu (11) była minimalna.

Podczas wyznaczania modelu rozmytego należy bardzo uważać na tzw. poziom walidacji poszczególnych modeli cząstkowych tzn. liczność zbioru danych, na podstawie którego został wyznaczony określony model lokalny θ_ij

(12)

Wartość ta nie może być zbyt mała, tzn. k_ij >> μ (μ liczba współczynników w wektorze modelu θ_ij), ponieważ w przeciwnym przypadku może zostać wyznaczony model θ_ij o błędnie określonej dynamice, np. może być niestabilny. Zwykle przyjmuje się relację k_ij ≥ (10-20)μ za bezpieczne spełnienie warunku dobrego określenia modelu θ_ij. Na poziom walidacji należy zwracać uwagę szczególnie przy większej liczbie sygnałów rozmywających, ponieważ w takim przypadku liczba modeli cząstkowych może być znaczna.

Modele parametryczne dynamiki protezy POLVAD

Parametryczny model dynamiki, odpowiadający schematowi oddziaływań z rys. 4, został oszacowany dla danych zarejestrowanych dla zadanej różnicy ciśnień ΔP = 100 mmHg. Parametry transmitancji (3) modelu (rząd n = 5, opóźnienia d_j należący do zbioru [0, 30 ms]) zostały optymalnie dobrane w sensie minimum błędu odtwarzania sygnału przepływu wylotowego F_wyl. Jakość tego modelu jest przedstawiona na rys. 6, gdzie zamieszczono pomierzony i wyznaczony na podstawie modelu wycinek przebiegu wypływu F_wyl.

Na podstawie powyższych przebiegów można uznać, że wyznaczony model zachowuje się całkiem poprawnie, a wyznaczony średni poziom błędów ε, określony jako ∑|ε| / Σ|F_wyl|, był równy 23,5 %. Obserwowany efekt wstecznych przepływów, przedstawiony na rys. 3, który niepokoił zespół badający zjawisko, został odtworzony przez model w całkiem satysfakcjonujący sposób, zarówno pod względem amplitudy zmian, jak czasu trwania.

Należy jednak pamiętać, że przedstawione wyniki modelowania (rys. 6), były uzyskane dla tej samej sekwencji danych, na podstawie których był wyznaczany model, tzn. dla zadanego spadku ciśnienia ΔP = 100 mmHg. W przypadku wykorzystania tego modelu dla innej sekwencji, pomierzonej dla zadanej różnicy ciśnień ΔP = 135 mmHg (rys. 7), model ten był wyraźnie błędny i wykazywał średni poziom błędów 157 % (!), tzn. model odtwarzał zmiany, ale w mniejszym zakresie i na całkiem innym poziomie.

Taki efekt oznaczał, że można „wpasować” model do określonej sekwencji danych, ale będzie on miał charakter bardzo lokalny i trzeba znaleźć wiele takich modeli w celu właściwego odtwarzania sygnału F_wyl.

W celu właściwego odtworzenia właściwości dynamicznych postanowiono wykorzystać modele rozmyte. Jako zmienne parametryzujące wykorzystano dwa sygnały: sygnał przepływu F_wyl oraz sygnał zadanej różnicy ciśnień ΔP. W przypadku sygnału przepływu F_wyl postanowiono przeprowadzić rozmycie uzależniające właściwości modelu od stanu: przepływ dodatni (właściwy wyrzut), przepływ wsteczny (niedomknięta zastawka) oraz faza domykania zastawki, co spowodowało ustawienie funkcji przynależności jak na rys. 8.

W przypadku drugiej zmiennej rozmywającej model podział na funkcje przynależności wyglądał podobnie (również trzy obszary), przy czym za niską różnicę ciśnień przyjęto wartości bliskie ΔP < 125 mmHg, a za wysokie ΔP powyżej 160 mmHg (rys. 9).

Przy takiej parametryzacji podczas wyznaczania modelu rozmytego należało określić 9 modeli cząstkowych. Pomiary zostały przeprowadzone dla różnych zadanych wartości ΔP = 100, 135, 150, 165, 175 mmHg. W każdym cyklu było 10 tys. pomiarów, a wyznaczane modele cząstkowe były określone na podstawie od ok. 300 do ponad 6000 pomiarów (12). W celu oszacowania wykorzystano 3 cykle: ΔP = 100, 150, 175 mmHg, a w celu weryfikacji dodatkowo ΔP = 135, 165 mmHg. Porównawcze wycinki przebiegów zostały przedstawione na rys. 10. Można zaobserwować wyraźną zgodność przebiegów mierzonych i modelowanych. Błędy modelowania były dla sekwencji ΔP = 135 mmHg, największe i wynosiły 23,3 %, a dla sekwencji ΔP = 175 mmHg tylko 10,4 %. Należy przy tym zwrócić uwagę, że występujące błędy są w większości efektem niedoskonałości techniki pomiarowej.

Model dobrze odtwarzał przebiegi zarówno dla danych, na podstawie których był wyznaczany (ΔP = 100, 150, 175 mmHg) jak i dla danych weryfikacyjnych (ΔP = 135, 165 mmHg) i można go uznać za dobre odtworzenie zachowania badanego procesu – zachowania dynamicznego komory POLVAD.

W artykule nie będą dyskutowane poszczególne transmitancje modeli cząstkowych, bo jest ich dużo: 9 modeli, każdy zawiera 5 transmitancji rzędu 4. W celu zobrazowania zmienności wyznaczonych modeli zostały przedstawione (rys. 11), odpowiedzi skokowe modeli przepływu F_wyl na zmianę ciśnienia P_zas w zależności od obszarów, dla których był wyznaczany model cząstkowy, przy czym na wykresach podano wskaźnik funkcji przynależności dla obydwu zmiennych rozmywających oraz wielkość wzmocnienia statycznego K modelu cząstkowego. Przykładowo, F_wyl = 2 i ΔP = 3 oznacza model cząstkowy wyznaczony dla: F_wyl zawartego w przedziale [-0,5; 0,5] oraz ΔP > 165.

Komentując wyniki (rys. 10) można uznać, że wyznaczony model nie obniża poziomu reakcji na pobudzenia, co dość często występuje w modelach określanych za pomocą estymatora LS, gdy są obecne zakłócenia pomiarowe, jak w rozważanym przypadku. Może być zatem wykorzystany do modelowania protezy POLVAD w różnych warunkach pracy i sprawdzania różnych, możliwych do rozważenia sposobów pobudzania. Takim krytycznym parametrem aplikacyjnym urządzeń tego typu jest szybkość zmian przepływu podczas wyrzutu (zbyt duży gradient zmian przepływu powoduje rozbijanie czerwonych krwinek).

Z tego powodu stosowane algorytmy sterowania oraz reakcje serwozaworów, sterujących zmianą ciśnienia P_a w części pneumatycznej protezy, są bardzo istotnym problemem. Wyznaczony model bardzo zgodnie (z pomiarami) odtwarza ten fragment mierzonych przebiegów i należy uznać, że będzie się bardzo dobrze nadawał do prowadzenia dalszych prac badawczych.

Przedstawione wybrane odpowiedzi skokowe (rys. 11) również wykazują wyraźne różnice. Należy zwrócić uwagę na zmienny współczynnik wzmocnienia K, który w przypadku ostatniego modelu cząstkowego ma wartość ujemną. Świadczy to o zmienności dynamiki w zależności od sygnałów rozmywających, jak również o potrzebie poszukiwania różnych modeli cząstkowych.

Podsumowanie

Przedstawiony w artykule sposób wykorzystania metody wyznaczania modeli rozmytych do określania właściwości dynamicznych procesów, dla których nie jest znany dokładny model działania, jest bardzo efektywnym podejściem i może w znaczący sposób przyspieszać badania, upraszczać technikę modelowania i przyspieszać etap wstępnej analizy działania całości układu, który zawiera zespół lub proces zastąpiony modelem rozmytym. Jednocześnie fakt, że model wynikowy pozwala zastępować działanie procesu w szerokim zakresie zmian punktu pracy procesu, umożliwia prowadzenie obliczeń optymalizacyjnych i szybką analizę działania alternatywnych rozwiązań układu lub diagnozowania procesu.

Podziękowania

Powyższa praca była wynikiem badań sponsorowanych przez narodowy program badawczy „Polskie Sztuczne Serce”, zgodnie z umową nr 12/WK/P02/0001/SPB-PSS/2009.

Bibliografia

1. Golnik A., Fajdek B., Janiszowski K.: Modelowanie układu krążenia z możliwością równoległego wspomagania pracy komór serca w pakiecie programowym PExSim. PAR 11/2010.
2. Janiszowski K.: Identyfikacja parametryczna układów dynamicznych w przykładach. Exit, 2002.
3. Janiszowski K. Wnuk P.: PexSim – a novel approach to the problem of the investigation of complex dynamic systems in an industrial environment. Problemy eksploatacji, 4/2006, s. 17–37.
4. Janiszowski K.: Modele procesów. „Modelowanie diagnostyka i sterowanie nadrzędne procesami, Implementacja w systemie Diaster”, Pr. zb. Red. J. Korbicz, J.M. Kościelny, WNT 2009.
5. Kozarski M., Ferrari G., Zielczyński K., Górczyńska K., Pałko K. J., Tokarz A., Darowski M.: A new hybrid electro-numerical model of the left ventricle. Computers in Biology and Medicine, vol.38, pp. 979–989, 2008.
6. Wnuk P.: Algorytmy optymalizacji struktury rozmytych modeli. Rozprawa doktorska, Warszawa 2004.

prof. dr hab. inż. Krzysztof Janiszowski

Działalność naukowa w zakresie identyfikacji, modelowania i sterowania procesów. Autor 5 podręczników, 35 artykułów, 84 referatów i ponad 40 projektów badawczych. Obszary badań aplikacyjnych: identyfikacja procesów przemysłowych (pieca obrotowego w cementowni, pieca do ciągłego wytopu stali, reaktorów chemicznych, pieca pirolitycznego, walczaka parowego, zespołu oczyszczalni itd.), synteza algorytmów sterowania do pneumatycznych układów pozycjonujących, hydraulicznych napędów wyporowych, zespołu napędu hybrydowego o 2 stopniach swobody, oraz pakiety oprogramowania do identyfikacji układów dynamicznych (IDCAD), i modelowania, symulacji działania i sterowania procesów przemysłowych (PExSim).