Simulink: równania różniczkowe

W niniejszym wpisie zakładam że czytelnikowi znane jest pojęcie „równanie różniczkowe” i całkowanie. Pierwotna wersja wpisu miała powstać przy użyciu XCOS wbudowanego w scilab. Niestety pojawiły się problemy. Uściślając za modelowanie prostego równania różniczkowego nie stanowiło problemu, aczkolwiek pojawiły się bardzo poważne problemy z stabilnością xcos. O ile sam scilab jest dość stabilny to moduł xcos już nie, wyłączanie się programu po symulacji czy w trakcie zdarzało się bardzo często na linuksie, a na windows już nie co rzadziej ale nadal było to normą. Niestety przyczyny mi nie są znane zważywszy na to że linuks korzystał z Openjre 6 a Windows oracle jre 7, zaryzykuję stwierdzenie że program jest stabilniejszy na wydaniu jre oracle ,aczkolwiek zauważmy oracle zapowiedział oparcie swoich wydań jre o openjre. Dodatkowo po za modelowaniu przykładu drugiego okazało się że wyniki nie są takie jak powinny być. Jeśli znajdzie się osoba która za modeluje przykład drugi w xcos i wyniki będą się zgadzały z wynikami z simulinka proszę o podesłanie pliku.

„...Mathematics is the language of nature...”

Powyższy cytat pochodzi z filmu PI i niesie w sobie niezwykle istotną informację. Zjawiska fizyczne mogą być opisane za pomocą równań(mniej lub bardziej dokładnie). Czytelniku jeśli nie wierzysz porównaj kształt huraganu z tzw. złotym podziałem.

Masa-Tłumik-Sprężyna

Zaczniemy od prostego przykładu układu II rzędu, który to w tego typu opracowaniach jest oklepany jak chabeta po westernie :) W poniższym dokumencie przykład 1 Link Przekształćmy równanie nieco(apostrofy to kolejne pochodne). F=mx''+cx'+kx F-siła m-masa c-stała tłumienia k-wsp. Sprężystości Nasuwa się pytanie jak przekształcić takie równanie aby można było je zapisać w simulinku. Otóż po jednej stronie musi zostać najwyższa pochodna a reszta po drugiej czyli: x''=1/m(F-cx'-kx) (1.1) Zapiszmy to teraz w simulinku:

Źródłem siły F jest blok Step czyli w chwili równej t=0, Pojawia się stały w czasie sygnał o wartości 1(po kliknięciu lpm, można zmienić wartość), następnie jest on przekazywany do sumatora, aby stworzyć odpowiedni sumator przeciągamy go z biblioteki(Math Operations) klikamy lpm po czym w list of signs wpisujemy |-+- ,a icon shape na rectangual. Przejdźmy do implementacji równania, w równaniu (1.1) mnoży się współczynniki w nawiasie przez odwrotność masy, w takim wypadku musimy skorzystać z bloku Gain, klikamy lpm na niego i wpisujemy po prostu 1/0.1 ,gdzie 0.1 to masa. Teraz każdy sygnał jaki do niego wejdzie(z sumatora) na wyjściu jest pomnożony przez ta wartość, następnie wyjście tego bloku musimy połączyć z Integratorem(będziemy to całkować), na wyjściu tego integratora pojawia się wtedy x' co musimy pomnożyć przez stałą c, czyli po raz kolejny Gain którego wejściem jest wyjście pierwszego integratora a wyjście wędruje do sumatora(nic innego jak sprzężenie zwrotne), dodatkowo z wyjścia pierwszego integratora ciągniemy do wejścia drugiego integratora na którego wyjściu pojawi się x czyli nasze rozwiązanie(konkretni położenie), z mnożeniem postępujemy tak jak wcześniej. Aby wykreślić rozwiązanie korzystamy z bloku Scope, po zakończeniu symulacji klikamy na niego.

A odpowiedź to:

Proponuje w celu polepszenia wyniku pobawić się solverem(ctr+e lub opcje konfiguracji).

Pokój

W niniejszym przykładzie zostanie pokazane krzywa temperatury nagrzewania się pokoju, zostanie także użyty regulator PI w celu polepszenia charakterystyki. Dlatego za nim pójdziesz dalej czytelniku zapoznaj się z regulatorami PID, opis znajduje się na Wikipedii a ponadto wiele uczelni udostępnia PDF'y wraz z opisem tego typu urządzeń.

Zacznijmy od implementacji równania, zauważmy że mamy tu układ równań. Stworzenie opisu w simulinku takich równań jest dość łatwe i podobne do poprzedniego. Używamy takich samych bloków jak poprzednio, zapewne zauważyłeś że wartości których szukamy są ze sobą połączone(Tg-temp. Grzejnika, Tp-temp. pomieszczenia) czyli mówiąc najoględniej wyjścia(czyli rozwiązania) są wejściami niektórych węzłów sumatorów. Zauważyłeś także zapewne że nie ma skąd pobrać wartości (Pg-moc grzania,tzew-temp. Zewnętrzna), dane te będą pochodziły z zewnątrz a nasza implementacja równań zostanie zamknięta w tzw. subsystemie czyli stworzymy tka jak by własny blok z 2 wejściami i wyjściami.

W bibliotece sources odnajdujemy blok In1(tylko jedno wyj.), potrzebne są dwa jeden z nich to zmienna Pg a drugi tzew. Teraz potrzebne są nam wyjścia Tp i Tg. W Sinks znajduje się blok Out1(tylko jedno wej.), potrzebne są dwa. Wyjście In łączymy do odpowiednich wejść sumatora a wyjścia do wyjść integrtorów. Teraz zaznaczamy wszystko i klikamy prawym przyciskiem i wybieramy create subsytem. Tada !

Mając pojedynczy blok opisujący równania Przystąpmy do stworzenia systemu sterowania. Potrzebny będzie blok Step symbolizujący zadaną temperaturę, sine wave(sinusoida) reprezentujący zmianę temperatury na zewnątrz(lpm i ustawiamy amplitudę na 5), Saturation który ograniczy sygnał(sygnał to moc grzejnika) ustawiamy wartość na 5000 a na koniec blok PI, możemy go zbudować sami(tak jak ma to zrobione) lub skorzystać z wbudowanego PID, aby otrzymać PI wartość D ustawiamy na zero. Skorzystamy także z bloku to workspace(wysłanie do matlaba) aby za pomocą polecenia plot wykreślić przebiegi(można też bezpośrednio w simulinku użyć scope i multipleksera), ustaw w save format bloku Array. Teraz łączymy wszystko jak na rysunku i symulujemy. Ustawienia symulacji to Type:fixed-step i ode1(Euler), 7000 sekund.

Celowo nie podaje nastaw PI. Spróbuj się pobawić czytelniku lub po prostu użyj jakieś metody doboru współczynników. Poniżej jest wykres w którym przedstawiam wyniki dla dwóch przykładów kiedy mamy do czynienia z regulatorem P tylko(nastawa P=5000) i regulatorem PI.