Vysoká škola báňská Technická univerzita Ostrava Fakulta strojní. Hypertextová podpora výuky v oblasti automatického řízení

Vysoká škola báňská Technická univerzita Ostrava Fakulta strojní

Hypertextová podpora výuky v oblasti automatického řízení

Vedoucí bakalářské práce: Ing. Renata Wagnerová Ph.D. Diplomant: Lukáš Hédl Ostrava: 2003

ANOTACE BAKALÁŘSKÉ PRÁCE HÉDL, L. Hypertextová podpora výuky v oblasti automatického řízení, Ostrava: katedra ATŘ – 352 VŠB-TU, 2003, 70 s. bakalářská práce, vedoucí Wagnerová, R. Bakalářská práce se zabývá tvorbou elektronické příručky pro podporu výuky v oblasti automatického řízení. Tato příručka obsahuje vybrané typové příklady z oblasti automatického řízení řešené pomocí programového systému MATLAB. Příručka obsahuje tři hlavní části. V první části je stručně popsán programový systém Matlab a jeho toolboxy. Další dvě části obsahují řešené příklady z oblasti analýzy a syntézy regulačních obvodů. Příručka je zpracována jako HTML dokument a bude dostupná na internetu pro potřeby studentů.

ANNOTATION OF BACHELOR THESIS HÉDL, L. The hypertext teaching support in the area of automatic control, Ostrava: Department of Control System and Instrumentation, Technical University of Ostrava, 2003, 70 p. Bachelor thesis, head: Wagnerova, R. The bachelor thesis is concerning with a creation of electronic manual for the teaching support in the field of automatic control. This manual contains selected model examples from the field of automatic control which are solved by the system MATLAB. The manual is divided into three parts. In the first part the programme system MATLAB and his toolboxes are described. Other two parts contain solved examples from the area of analysis and synthesis of the control circuits. The manual is processed as a HTML document and it will be accessible for the use of students on the internet.

Obsah

1

Obsah Seznam použitých symbolů................................................................................................. 2 1

Úvod .............................................................................................................................. 4

2

Popis systému Matlab a jeho toolboxů ...................................................................... 5 2.1 Matlab.................................................................................................................... 5 2.2 Simulink ................................................................................................................ 6 2.3 Control System Toolbox........................................................................................ 7 2.4 System Identification Toolbox .............................................................................. 7 2.5 Optimization Toolbox ........................................................................................... 7 2.6 Signal Processing Toolbox .................................................................................... 8

3

Analýza systémů .......................................................................................................... 9 3.1 Zadání matematického modelu systému................................................................ 9 3.1.1 Diferenciální rovnice..................................................................................... 9 3.1.2 Přenos systému ............................................................................................ 18 3.1.3 Zadání přenosu systému pomocí pólů a nul ................................................ 20 3.1.4 Zadání stavového popisu ............................................................................. 22 3.1.5 Převody mezi jednotlivými způsoby popisu systému.................................... 24 3.2 Propojení systémů ............................................................................................... 28 3.2.1 Paralelní zapojení ....................................................................................... 28 3.2.2 Sériové zapojení........................................................................................... 29 3.2.3 Zpětnovazební zapojení ............................................................................... 29 3.3 Vykreslení základních charakteristik systému .................................................... 31 3.3.1 Přechodová charakteristika......................................................................... 31 3.3.2 Impulsní charakteristika .............................................................................. 34 3.3.3 Kmitočtové charakteristiky .......................................................................... 36 3.4 Stabilita systému.................................................................................................. 39 3.4.1 Stabilní systém ............................................................................................. 40 3.4.2 Systém na mezi stability............................................................................... 43 3.4.3 Nestabilní systém ......................................................................................... 46

4

Syntéza regulačních obvodů ..................................................................................... 49 4.1 Návrh regulátoru metodou Ziegler – Nichols...................................................... 49 4.2 Návrh regulátoru metodou požadovaného modelu.............................................. 52 4.3 Metody výpočtu koeficientů diskrétních regulátorů............................................ 58

5

Závěr ........................................................................................................................... 63

6

Literatura ................................................................................................................... 64

Seznam použitých symbolů

Seznam použitých symbolů A

stavová matice systému typu (n,n)

A(ω)

modul, amplituda kmitočtového přenosu

Ao

amplituda otevřeného regulačního obvodu

ai

koeficient mnohočlenu ve jmenovateli přenosu, koeficient levé strany lineární diferenciální rovnice

b bj

stavová matice řízení typu (n,1) koeficient mnohočlenu v čitateli přenosu, koeficient pravé strany lineární diferenciální rovnice

ci

koeficienty homogenní diferenciální rovnice

c E

výstupní matice systému typu (n,1)

f0

frekvence [s-1]

g

impulsní funkce

G(jω)

kmitočtový přenos

G(s)

obrazový přenos (L – přenos)

G(z)

obrazový přenos diskrétního systému (Z – přenos)

Go

přenos otevřeného regulačního obvodu

GMČ

přenos měřicího členu

GP

přenos, přes který působí na regulační obvod poruchová veličina

GR

přenos regulátoru

GS

přenos regulované soustavy

Gv

přenos poruchy

Gve

odchylkový přenos poruchy

Gw

přenos řízení

Gwe

odchylkový přenos řízení

h

přechodová funkce

Im{G(jω)}

imaginární část kmitočtového přenosu

j = −1

obraz regulační odchylky

imaginární jednotka

k1

koeficient přenosu

L

operátor přímé Laplaceovy transformace

L(ω)

logaritmická amplituda kmitočtového přenosu [dB]

-1

L

operátor zpětné Laplaceovy transformace

m

stupeň mnohočlenu v čitateli kmitočtového přenosu

n

stupeň mnohočlenu ve jmenovateli přenosu

No

charakteristický mnohočlen otevřeného regulačního obvodu

2

Seznam použitých symbolů P(ω)

reálná část kmitočtového přenosu

Q(ω)

imaginární část kmitočtového přenosu

Re{G(jω)}

reálná část kmitočtového přenosu

s

komplexní proměnná Laplaceovy transformace

t

spojitý čas [s]

T, Ts

vzorkovací perioda

Td

dopravní zpoždění u spojitých systémů (členů)

TD

derivační časová konstanta

TI

integrační časová konstanta

U

obraz vstupní veličiny, obraz řídicí (akční) veličiny

u

vstupní veličina, řídicí (akční) veličina

V

obraz poruchové veličiny

v

poruchová veličina

W

obraz žádané veličiny

w

žádaná veličina, nezávisle proměnná u bilineární transformace

x Y

vektor stavových veličin (stav) dimenze n

y

regulovaná (výstupní) veličina

z

nezávisle proměnná Z transformace

obraz regulované (výstupní) veličiny

Z -1

operátor přímé transformace Z

Z

operátor zpětné transformace Z

ω

reálný úhlový kmitočet [s-1]

ξ

koeficient poměrného tlumení

δ

Diracův jednotkový impuls

η

Heavisideův jednotkový skok

ϕ(ω)

fáze kmitočtového přenosu [rad]

3

Úvod

4

1 Úvod V dnešní době existuje množství programů usnadňujících matematické výpočty z oblasti automatického řízení. Patří mezi ně i programový systém Matlab. Právě tento systém je jedním z programových produktů využívaných na katedře ATŘ. Síťová verze tohoto programu je dostupná studentům na školních počítačových učebnách. Proto vznikla potřeba elektronické příručky, která by studenty seznámila s možnostmi tohoto systému. Cílem této bakalářské práce je zpracování elektronické příručky pro podporu výuky automatického řízení. Tato příručka obsahuje řešení typových příkladů z oblasti automatizace pomocí programového prostředí Matlab. Příručka je rozdělena na tři hlavní části. První kapitola obsahuje seznámení s programem Matlab a jeho toolboxy Simulink, Control System, System Identification, Signal Processing a Optimization, také popisuje možnosti jejich využití. Druhá kapitola je zaměřena na oblast analýzy systémů, možnosti zadávání matematických modelů a převodů mezi nimi, vykreslení základních charakteristik systémů a kontrolu stability. Třetí část obsahuje příklady z oblasti syntézy regulačních obvodů jsou tu uvedeny postupy pro návrh spojitých a diskrétních regulátorů v prostředí Matlab. Elektronická podoba příručky je vytvořena jako HTML dokument a je umístěna na internetových stránkách katedry ATŘ spolu s ostatními již vytvořenými učebnicemi.

Popis systému Matlab a jeho toolboxů

5

2 Popis systému Matlab a jeho toolboxů 2.1 Matlab Matlab je interaktivní systém pro vědecké a technické výpočty založený na maticovém kalkulu. Umožňuje řešit velkou oblast numerických problémů, aniž byste museli programovat vlastní program. Název Matlab vznikl zkrácením z MATrix LABoratory. Původně Matlab vznikl jako interaktivní nadstavba pro usnadnění práce s knihovnami LINPACK a EISPACK pro práci s maticemi. Základním typem byly matice, které na rozdíl od většiny jiných systémů a jazyků nevyžadovaly nastavování dimenzí. Současný Matlab je mnohem více než jen nadstavbou maticové knihovny.Systém obsahuje vlastní interpretr jazyku Matlab, ve kterém lze připravit jak dávkové soubory, tak definovat i nové funkce. Tyto funkce mohou být interpretovány buď přímo z textové podoby souborů nazývaných m-file nebo z předzpracované podoby p-file. Systém dále umožňuje přidávat moduly (soubory mex-file) zkompilované do strojového kódu procesoru. Jazykem zdrojových souborů může být jazyk C, C++, Fortran nebo po použití mcc (Matlab to C Compiler) i funkce uložená v m-file. V prostředích, kde je možný grafický výstup, je k dispozici velmi silná podpora pro tvorbu uživatelského prostředí a vizualizaci dat. Základní funkce umožňují vizualizaci 2D a 3D dat v grafech s množstvím volitelných parametrů, které lze po doplnění legendou a popisy os snadno vytisknout. Systém umožňuje i tvorbu vlastních dialogů a oken s kombinacemi grafů, tlačítek, listů a dalších objektů, kterým mohou být přiřazeny volané funkce v m-file. Dobře je obslouženo i polohovací zařízení - myš. Pro tvorbu samostatně spustitelných aplikací je k dispozici Matlab library. Asi nejdůležitější částí instalace Matlab jsou ''knihovny'' funkcí (ve skutečnosti adresáře s m a mex soubory), které jsou nazývány toolboxy. Toolboxy obsahují vždy uceleným způsobem včetně dokumentace a příkladů zpracovaný určitý obor numerické matematiky, analytické matematiky, statistiky, systémového přístupu k regulacím a další obory, ve kterých nachází Matlab uplatnění. Po výčtu základních vlastností systému Matlab je lépe možné si představit uplatnění systému v typických oblastech pro jeho použití tak, jak je udává návod k systému Matlab: • Matematické výpočty • Vývoj algoritmů • Modelování a simulace • Analýza dat a vizualizace • Vědecká a inženýrská grafika • Vývoj aplikací včetně uživatelského interface

Popis systému Matlab a jeho toolboxů

6

2.2 Simulink Simulink postupně přerostl z knihovny funkcí určené k simulaci jednoduchých lineárních spojitých a diskrétních systémů v samostatný subsystém s dokonalým uživatelským rozhraním. Základem toolboxu Simulink jsou bloky, které reprezentují elementární dynamické systémy. Propojením signálových vstupů a výstupů těchto bloků vznikají modely složitějších systémů. Libovolnou skupinu bloků lze uzavřít do subsystému a určit externí vstupy a výstupy této skupiny. Dále lze pracovat s takovouto skupinou jako se základním blokem. Je-li potřeba zastínit proměnné parametry bloků uvnitř skupiny, lze uzavřenou skupinu zamaskovat a doplnit informacemi, které vytvoří při modifikaci parametrů bloku dotazový dialog a postarají se o přepočítání a přenesení zadaných parametrů dovnitř do zamaskované skupiny. Pro zamaskovanou skupinu lze také vytvořit grafickou reprezentaci skupiny, která se může být i závislá na nastavených parametrech. K výpočtům parametrů lze užít všech forem výrazů a volání funkcí, které Matlab umožňuje. Interaktivní způsob tvorby a simulace modelů se spouští z příkazové řádky systému Matlab příkazem simulink. Po spuštění je vytvořeno okno pro tvorbu nového modelu a okno obsahující základní nabídku otvírání knihoven zdrojů signálů, základních spojitých, diskrétních a nelineárních bloků a bloků pro zobrazování a ukládání signálů. Pod touto interaktivní obálkou se skrývá systém velmi podobný grafickému subsystému s obdobnými funkcemi simget a simset. Další vrstva funkcí umožňuje již komfortnější neinteraktivní tvorbu modelů systémů. Pro obvyklého uživatele však není nutné o implementaci a programování modelů přemýšlet. Základní a již dále nedělitelné jsou pouze bloky obsahující takzvané s-funkce. Jedná se o zabudované funkce, mex-soubory a nebo o obvyklé interpretované funkce uložené v m-souborech. Tyto funkce mají předepsané parametry a chování pro různé druhy volání. Pro spojité systémy v jednoduchosti informují o okamžitých hodnotách derivací a výstupů bloku pro diskrétní bloky o příštích hodnotách výstupů při zadaných vstupních hodnotách. Jednoduché nelineární bloky lze také vytvářet vložením libovolného přípustného výrazu v jazyku Matlab do k tomu určeného bloku. Simulink je schopen simulovat smíšené systémy obsahující spojité části, diskrétní části i s různými periodami vzorkování a s posunutými okamžiky vzorkování. Je schopen simulovat i nelineární bloky a aproximovat chování systémů obsahujících algebraické smyčky, na které ovšem před simulací upozorňuje. Dynamické vlastnosti lineárních částí lze popisovat komplexními přenosy, maticemi systémů nebo přímo použít bloky reprezentující přímo sčítání, integraci, diferenci, násobení konstantou a další elementární operace. V knihovně nelineárních bloků jsou předdefinovány paměťové bloky, přepínače, releové charakteristiky, násobení a dělení signálů, zdroje hodinových impulsů a mnoho dalších. Většinu systémů, které mohou být řešeny v diskrétním čase (lze je diskretizovat) je možno přes RTW (Real Time Workshop) převést přímo do zdrojového kódu v jazyku C, který lze po doplnění zdrojovými texty pro uživatelem definované s-funkce přímo kompilovat do strojových jazyků počítačů určených k řízení. Pak vytvoření i složitého regulačního systému může vypadat tak, že návrhář sestaví model navrhovaného regulátoru z předdefinovaných bloků Simulink a po stisknutí jediného tlačítka (je-li již správně zvolena konfigurace) dojde ke kompletnímu překladu a přenosu výsledného kódu do řídící jednotky.

Popis systému Matlab a jeho toolboxů

7

2.3 Control System Toolbox Control System Toolbox je aplikační knihovna, která rozšiřuje systém Matlab o nástroje pro řídicí techniku a teorii systémů. Funkce z oblasti analýzy a návrhu řídicích systémů využívají jak klasické přechodové charakteristiky, tak i popisy systémů ve stavovém prostoru. Novinkou je zavedení lineárních časově invariantních objektů (LTI), což jsou struktury popisující jednorozměrové i mnoharozměrové lineární systémy. Do LTI lze kromě popisu struktury systému uložit i mnoho dalších vlastností, jako je vzorkovací frekvence, dopravní zpoždění, pojmenování vstupních a výstupních signálů a další uživatelská data. Tyto informace lze samozřejmě editovat a ukládat v kterémkoli časovém okamžiku a tak přehledně dokumentovat jednotlivé stavy systému během experimentu. Operace s LTI jsou podobné maticovým operacím (sčítání, násobení, …). LTI umožňuje uživateli pracovat s přenosy systémů, se stavovým prostorem i s popisy pomocí pólů a nul systému. Vestavěný grafický LTI Viewer poskytuje nástroje na analýzu odezvy systému, jako jsou přechodová charakteristika, kmitočtové charakteristiky v logaritmických souřadnicích i komplexní rovině, zobrazení pólů a nul a další. Pouhým klepnutím na tlačítko myši je možno přecházet z časové do kmitočtové oblasti, volit množinu pozorovaných vstupů a výstupů, pozorovat pouze podstatné části charakteristik pomocí funkce "zoom" a podobně.

2.4 System Identification Toolbox System Identification Toolbox je určen pro vytváření matematických modelů systémů z naměřených dat. Poskytuje nástroje pro vytvoření matematických modelů dynamických systémů založené na sledování vstupních a výstupních dat. Pro práci využívá grafické uživatelské rozhraní, které usnadňuje práci při organizaci dat a modelů. Nástroje pro identifikaci poskytované v tomto toolboxu jsou určené pro užití v oblasti od návrhu regulátoru a zpracování signálu až k časové a vibrační analýze.

2.5 Optimization Toolbox Optimization toolbox je určen pro minimalizaci a maximalizaci funkcí. Obsahuje funkce určené pro minimalizaci (nebo maximalizaci) obecných nelineárních funkcí. Také obsahuje funkce pro nelineární vyrovnání řešení a funkce pro řešení problému nejmenších čtverců.

Popis systému Matlab a jeho toolboxů

8

2.6 Signal Processing Toolbox Tento nástroj je určen ke zpracování signálů. Podporuje široké pole operací ke zpracování signálů od generování časových průběhů signálu po návrh filtrů a jejich implementaci, parametrické modelování a spektrální analýzu. Toolbox nabízí dvě kategorie nástrojů: •

Funkce pro zpracování signálů

•

Grafické, interaktivní nástroje

Funkce pro zpracování signálů jsou volány přímo z příkazové řádky v prostředí Matlab. Tyto funkce jsou většinou naprogramovány jako m-file a po zkopírování a přejmenování je možné je modifikovat. Grafické uživatelské rozhraní (GUI) nabízí prostředí pro návrh filtrů, analýzu signálů a jejich implementaci, nástroje pro prohlížení průběhu signálů a jejich editaci. Grafické prostředí umožňuje prácí s myší a grafickou editaci signálů, signály je možné přehrát na zvukovém zařízení počítače, a mnoho dalších.

Analýza systémů

9

3 Analýza systémů 3.1 Zadání matematického modelu systému Dynamické vlastnosti systému lze popsat různými způsoby, které lze v podstatě rozdělit do dvou skupin a to na vnitřní a vnější popis systému. Vnější popis systému je vyjádření dynamických vlastností systému závislostí mezi jeho vstupem a výstupem. Při vnějším popisu systému považujeme systém za černou skříňku se vstupem a výstupem. Vnějších popisů lineárního systému je několik. Všechny systém jednoznačně charakterizují a povětšině lze z jednoho druhu popisu odvodit další. Pod pojmem popis lineárního regulačního systému si nejčastěji představíme popis řízeného systému (regulované soustavy). Tomu napovídá i používané označení u(t) pro vstupní veličinu a y(t) pro výstupní veličinu. Ale můžeme tím myslet i popis prvků systému. Rovněž tak můžeme stejným způsobem popisovat i řidicí systém, kde je opačné značení vstupní a výstupní veličiny. Vnější popis lineárního spojitého systému s jednou vstupní a jednou výstupní veličinou může být: • diferenciální rovnice • přenos a poloha pólů a nul přenosu systému • přechodová funkce • impulsní funkce • kmitočtový přenos Jediným vnitřním popisem systému je jeho popis ve stavovém prostoru tzv. stavový popis. Vnitřní popis systému je vyjádření dynamických vlastností systému závislostí mezi jeho vstupem, stavem a výstupem. Vnitřní popis systému je popis jeho stavově přechodné struktury Je nejdokonalejším popisem systému. Vyjadřujeme jej pomocí stavových rovnic. [Švarc, I. 1992]

3.1.1 Diferenciální rovnice Popis systému v časové oblasti. Lineární spojitý systém S s jednou vstupní a jednou výstupní veličinou podle obr. 3.1 má obecnou závislost mezi vstupem a výstupem popsán rovnicí ve tvaru d n y(t) d n−1y(t) dy(t) d mu(t) d m−1u(t) du(t) an + an−1 n−1 + ...+ a1 + a0 y(t) = bm m + bm−1 m−1 + ...+ b1 + b0u(t), (3.1) n dt dt dt dt dt dt Kde jsou: ai, bi - koeficienty levé a pravé strany diferenciální rovnice, u(t) - vstupní veličina, y(t) - výstupní veličina.

Analýza systémů

10 u(t)

y(t)

S

obr. 3.1 Blokové schéma systému

Diferenciální rovnici systému získáme obvykle tak, že uvedeme fyzikální vztahy systému a eliminujeme všechny proměnné mimo vstupní a výstupní veličinu. Např. v mechanických systémech obvykle vystačíme s dynamickou rovnováhou sil [Švarc, I. 1992]. Matlab umožňuje diferenciální rovnice řešit numericky. Pro zadání diferenciální rovnice do programu Matlab máme dvě možnosti. Rovnici můžeme zadat pomocí toolboxu Simulink, kde sestavíme její blokové schéma. Druhou možností je naprogramovat si pomocí m-file editoru vlastní funkci, do které zadáme parametry rovnice.

Příklady: Příklad 1. – Nádrž s volným odtokem kapaliny Provedeme analytickou identifikaci nádrže s volným odtokem kapaliny obr. 3.2 za předpokladu, že h(t) je výška hladiny [m], P je obsah příčného průřezu hladiny [m2], A je obsah příčného průřezu odtoku [m2], q1 je objemový přítok [m3s-1], q2 je objemový odtok [m3s-1] a ϕ rychlostní součinitel. Pro q2 platí: q 2 = Aϕ 2 gh(t )

(3.2)

Máme zadány následující hodnoty: P = 2 m2 q1 = 0,1 m3s-1 A = 0,1 m2 ϕ = 0,9 h(0) = 0,1 m q1(t)

P

h(t)

q2(t) obr. 3.2 Nádrž s volným odtokem

Pro elementární přírůstek objemu kapaliny v nádrži Pdh(t) za elementární časový přírůstek dt platí rovnice Pdh(t) = q1(t)dt – q2(t)dt, Resp. Po uvažování vztahu (3.2) Pdh(t) = q1(t)dt - Aϕ 2 gh(t ) .

Analýza systémů

11

Po úpravě dostaneme nelineární nehomogenní diferenciální rovnici 1. řádu Ph′(t ) = q1 (t ) − Aϕ 2 gh(t )

⇒

[

] P1

h′(t ) = q1 (t ) − Aϕ 2 gh(t ) ⋅

(3.3)

a) Řešení pomocí toolboxu Simulink Na obr. 3.3 je vidět blokové schéma diferenciální rovnice vytvořené v toolboxu Simulink, při sestavení blokového schéma vycházíme z rovnice (3.3).

obr. 3.3 Řešení diferenciální rovnice nádrže v programu Simulink

Popis jednotlivých bloků a jejich nastavení: Step – skoková změna pro nás představuje přítok nádrže q1. Zde nastavíme tyto parametry: Step time – zde uvádíme v jakém čase od začátku simulace dojde k zapnutí skokové změny čas uvádíme v sekundách. Hodnotu nastavíme na nulu. Initial value – počáteční hodnota. Hodnotu nastavíme na nulu. Final value – výsledná hodnota skoku kterou nastavíme na hodnotu 0,1. Gain – zesílení tento blok vynásobí vstupní signál požadovanou hodnotou. V tomto řešení je zesílení použito dvakrát: Gain – Zde je zadána hodnota 1 / P ve tvaru 1 / 2. Gain2 – V tomto zesílení je zadán vztah Aϕ 2 g ve tvaru 0.1*0.9*sqrt(2*9.81). Integrátor – Integrační člen v tomto členu zadáváme počáteční podmínku 0,1. Math Function – tento blok umožňuje provedené několika matematických funkci. Pro nás je důležitá funkce odmocnina – sqrt. Scope – nástroj na zobrazení výsledného signálu v závislosti na čase. Graf otevřeme poklepáním na tento blok. b) Řešení pomocí m-file Zadání funkce: function hp = difrov (t, h) % deklarace funkce a načtení hodnot: t-doba simulace, h-počáteční podmínka q1 = 0.1; %zadání konstant A = 0.1; %zadání konstant fi = 0.9; %zadání konstant g = 9.81; %zadání konstant hp = (q1 - A*fi*sqrt(2*g*h))/2; %zadání diferencialní rovnice

Analýza systémů

12

Název uloženého souboru se musí shodovat s názvem funkce. V tomto případě se soubor bude jmenovat difrov.m. Při volbě názvu funkce nesmíme zapomenout, že naše funkce nesmí mít stejný název jako vestavěné funkce programu Matlab. Abychom mohli vytvořenou funkci v podobě m-file v spustit musíme v programu nastavit cestu do adresáře kde je náš soubor uložen. Matlab standardně nabízí složku MATLABR11\work. Nastavení cesty ke složce můžeme změnit pomocí nástroje Path Browser (File -SetPath…), pomocí tohoto nástroje také můžeme cestu k našemu adresáři uložit pro příští spuštění. Výpočet: [t,y]=ode45('difrov', [0:0.1:30], 0.1); %volaní funkce difrov a výpočet

Funkce ode45 je numerická funkce řešení diferenciálních rovnic. Parametry [t, y] určují jména proměnných kde bude uložen výsledek. Parametr 'difrov' volá naši vytvořenou funkci. Hodnoty [0:0.1:30] určují krok a dobu řešení. Poslední uvedená hodnota 0,1 je počáteční podmínka y0. vykreslení: plot (t, y); %vykreslení grafu grid; %zapnutí zobrazení mřížky title ('Výška hladiny v nádrži'); %titulek grafu xlabel ('doba simulace'); %popis osy x ylabel ('výška hladiny'); %popis osy y

Výsledek:

obr. 3.4 Vykreslení přechodové charakteristiky modelu nádrže

Analýza systémů

13

Příklad 2. – RC člen Provedeme analytickou identifikaci elektrického obvodu se dvěma prvky kondenzátorem C a odporem R a vstupním napětím u1. Zapojení obvodu je patrné z obr. 3.5. Máme zadány hodnoty: R = 10 kΩ C = 10 µF u1 = 220V R

U1

C

U2

obr. 3.5 Zapojení RC členu

Při řešení elektrických obvodů vycházíme z Ohmova zákona a Kirchhoffových zákonů. Pro úbytek napětí na odporu R platí vztah u (t ) = Ri (t ) ⇒ U ( s ) = RI ( s ) ,

(3.4)

pro úbytek napětí na kondenzátoru C platí: t

u (t ) =

1 i (τ )dτ C ∫0

⇒ U ( s) =

1 I (s) . Cs

Poměr výstupního napětí na vstupním je možné popsat vztahem U 2 (s) Z 2 (s) = , U 1 ( s ) Z1 ( s ) + Z 2 ( s )

(3.5)

(3.6)

kde Z1(s) představuje příčnou větev obvodu (kondenzátor C) a Z2(s) horní větev obvodu (odpor R). Dosazením vztahů (3.4) a (3.5) do rovnice (3.6) získáme 1 1 I ( s) U 2 (s) 1 Cs = = Cs = . RCs + 1 RCs + 1 1  U1 ( s)  R +  I (s) Cs Cs   Výsledkem je rovnice vyjadřující závislost výstupního napětí u2 na vstupním napětí u1 ve tvaru U 2 (s ) 1 = , U1 (s ) RCs + 1 tuto rovnici upravíme na následující tvar RCsU 2 ( s ) + U 2 ( s ) = U 1 ( s ) .

Analýza systémů

14

Pomocí zpětné Laplaceovy transformace rovnici převedeme z oblasti komplexní proměnné s do časové oblasti a vyjádříme nejvyšší derivaci: RCu& 2 (t ) + u 2 (t ) = u1 (t ) u& (t ) = [u1 (t ) − u 2 (t )]

1 RC a) Řešení pomocí toolboxu Simulink

(3.7)

Na obr. 3.6 je vidět blokové schéma diferenciální rovnice vytvořené v toolboxu Simulink při sestavení blokového schéma vycházíme z rovnice (3.7).

obr. 3.6 Řešení diferenciální rovnice RC členu v programu Simulink

Popis jednotlivých bloků a jejich nastavení: Step – zde nastavíme vstupní napětí u1: Step time = 0, Final value = 220, Gain – do zesílení nastavíme hodnotu 1/RC, Integrátor – Integrační člen v tomto členu zadáváme počáteční podmínku 0, Scope – nástroj na zobrazení výsledného signálu v závislosti na čase. Graf otevřeme poklepáním na tento blok. b) Řešení pomocí m-file Zadání funkce: function u2 = difrov2 (t, u) % deklarace funkce a na čtení hodnot: tdoba simulace, u-počáteční podmínka u2 u1=220; %zadání konstant R=10000; %zadání konstant C=0.01; %zadání konstant RC=R*C; %zadání konstant u2= (u1 - u)/RC; %diferencialní rovnice

Výpočet: [t,u]=ode45('difrov2', [0:0.1:800], [0]); % volání funkce difrov2

Funkce ode45 je numerická funkce řešení diferenciálních rovnic. Parametry [t, u] určují jména proměnných kde bude uložen výsledek. Parametr 'difrov2' volá naši vytvořenou funkci. Hodnoty [0:0.1:800] určují krok a dobu řešení. Poslední uvedená hodnota 0 je počáteční podmínka. Vykreslení: plot (t, u); %vykreslení grafu grid; %zobrazení mřížky title ('průběh napětí u2'); %titulek grafu xlabel ('doba simulace'); %popis osy x ylabel ('napětí') % popis osy y

Analýza systémů

15

Výsledek:

obr. 3.7 Vykreslení přechodové charakteristiky modelu RC členu

Příklad 3. – Hmota na pružině a tlumiči Provedeme analytickou identifikaci systému s hmotou na pružině a tlumiči (obr. 3.8) za předpokladu, že k – tuhost pružiny [kg.s-2], b – odpor tlumiče [kg.s-1], m – hmotnost [kg] a F – působící vstupní síla [N]. Sledované výstupní veličiny jsou y′(t) – průběh rychlosti hmoty [m.s-1] a y(t) – průběh výchylky polohy hmoty [m]. Máme zadány hodnoty: F = 50N m = 2kg k = 5kg.s-2 b = 0,5kg.s-1 y0 = 5m y’0 = 0m.s-1

k

b

f(t)

m

y(t), y’(t)

obr. 3.8 Hmota na pružině a tlumiči

Analýza systémů

16

Při řešení vycházíme z pohybové rovnice m

ma = ∑ Fi , i =1

pro náš příklad má pohybová rovnice následující tvar my′′(t ) = F (t ) − Fk (t ) − Fb (t ) ⋅

(3.8)

Jednotlivé síly spočítáme podle následujících vztahů Fk = ky ; Fb = ky ′ . Po dosazení do rovnice (3.8) získáme diferenciální rovnici druhého řádu my ′′(t ) = f (t ) − ky (t ) − by ′(t ) , po úpravě získáme tvar y ′′(t ) = [ f (t ) − ky (t ) − by ′(t )]⋅

1 ⋅ m

(3.9)

a) Řešení pomocí toolboxu Simulink Na obr. 3.9 je vidět blokové schéma diferenciální rovnice vytvořené v toolboxu Simulink při sestavení blokového schéma vycházíme z rovnice (3.9).

obr. 3.9 Řešení diferenciální rovnice hmoty na pružině a tlumiči v programu Simulink

Popis jednotlivých bloků a jejich nastavení: Step – zde natavíme vstupní sílu F: Step time = 0, Final value = 50. Gain – hodnoty jednotlivých zesílení jsou patrné z obr. 3.9. Integrátor – Integrátor nastavíme na hodnotu 0, Integrátor1 na hodnotu 5 b) Řešení pomocí m-file Protože numericky je možné řešit pouze diferenciální rovnice prvního řádu, musíme nejdříve výchozí diferenciální rovnici druhého řádu (3.9) převést na soustavu diferenciálních rovnic prvního řádu. Nejdříve zavedeme substituci y1 = y; y2 = y’,

Analýza systémů

17

po zavedení substituce získáme soustavu dvou diferenciálních rovnic prvního řádu y1′ = y 2 y ′2 = [ f (t ) − ky1 − by 2 ] ⋅

1 m

(3.10)

Ze soustavy rovnic (3.10) získáme vektorový popis y´ = f - K . y  0  0 y´ =  f (t )  −  k  m   m

− 1  y  b ⋅ 1 y m   2 

(3.11)

Zadání funkce: function yp = difrov1 (t, y) % deklarace funkce a na čtení hodnot: tdoba simulace, y-počáteční podmínky f = [0;25]; % zadání konstant - Vektoru f k = [0 -1;2.5 0.25]; %zadání konstant - Vektoru K (hodnoty k, b) yp = f-k*y; %zadání diferenciální rovnice podle vztahu

Při zadávání konstant a diferenciální rovnice vycházíme ze vztahu (3.11) proto f a k zadáváme jako matice. Výsledná hodnota bude také matice yp, kde první sloupec bude obsahovat hodnoty y (průběh výchylky hmoty) a druhý sloupec hodnoty y’ (průběh rychlosti hmoty). Název uloženého souboru se musí shodovat s názvem funkce. V tomto případě se soubor bude jmenovat difrov1.m. Výpočet [t,yp]=ode45('difrov1', [0:0.1:50], [5;0]);

Funkce ode45 je numerická funkce řešení diferenciálních rovnic. Parametry [t, yp] určují jména proměnných, kde bude uložen výsledek. V proměnné yp bude uložen jako matice. Parametr 'difrov1' volá naši vytvořenou funkci. Hodnoty [0:0.1:30] určují krok a dobu řešení. Poslední hodnota [5;0] je matice počátečních podmínek. Vykreslení: plot (t, yp) %vykresleni grafu grid; %zobrazeni mřížky title ('průběh výchylky a rychlosti hmoty'); %titulek grafu xlabel ('doba simulace'); %popis osy x ylabel ('výchylka, rychlost'); %popis osy y

Analýza systémů

18

Výsledek:

y

y´

obr. 3.10 Vykreslení přechodové charakteristiky modelu hmoty na pružině a tlumiči

3.1.2 Přenos systému Popis systému v oblasti komplexní proměnné s. Přenos systému definujeme jako poměr Laplaceova obrazu výstupní veličiny k Laplaceovu obrazu vstupní veličiny při nulových počátečních podmínkách. Značíme ho symbolem G(s) a za předpokladu, že L{u(t)}=U(s) a L{y(t)}=Y(s) můžeme jeho definici psát ve tvaru G (s) =

L{u(t)} Y(s) = . L{y(t)} U(s)

(3.12)

Máme-li systém, popsaný diferenciální rovnicí (3.1), provedeme transformaci této rovnice za použití věty o linearitě a věty o obrazu n-té derivace v Laplaceově transformaci. Jestliže podle definice přenosu položíme všechny počáteční podmínky rovny nule, dostaneme po vytknutí na levé i pravé straně vztah (a n s n + ... + a1s + a0 )Y ( s ) = (b m s m + ... + b1s + b0 )U ( s ) odtud podle definice (3.12) vyplyne vzorec pro výpočet přenosu z koeficientů známé diferenciální rovnice G (s) =

Y(s) b m s m + ... + b1s + b0 = . U(s) a n s n + ... + a1s + a0

(3.13)

Tento vzorec má oproti definici v čitateli polynom, vytvořený ze vstupní strany diferenciální rovnice a ve jmenovateli polynom z výstupní strany diferenciální rovnice. [Švarc, I. 1992] Systém Matlab nám diky toolboxu Control System poskytuje několik možností zadání přenosu.

Analýza systémů

19

Příklady: Zadání spojitého přenosu za pomocí funkce tf - create a transfer function model Základní syntaxe příkazu: SYS = TF(NUM,DEN) – spojitý systém zadáváme čitatele a jmenovatele sys = tf(NUM, DEN,'OutputDelay',Td) – pokud potřebujeme zadat dopravní zpoždění přidáme parametr OutputDelay a hodnotu zpoždění Td. s = TF('s') – vytvoření komplexní proměnné s Máme systém zadaný přenosem: G (s) =

1,2 s + 0,5 7,1s + 2 s 2 + 1,5s + 0,1 3

Pro zadání přenosu použijeme příkaz: sys = tf([1.2 0.5], [7.1 2 1.5 0.1])

Do proměnné sys se uložil přenos v tomto tvaru: Transfer function: 1.2 s + 0.5 ----------------------------7.1 s^3 + 2 s^2 + 1.5 s + 0.1

Při zadávání přenosu můžeme nejdřív zadat do vlastních proměnných koeficienty čitatele a jmenovatele a teprve potom pomocí těchto proměnných zadat přenos: cit = [1.2 0.5]; jmen = [7.1 2 1.5 0.1]; sys = tf(cit, jmen) Transfer function: 1.2 s + 0.5 ----------------------------7.1 s^3 + 2 s^2 + 1.5 s + 0.1

Přenos můžeme zadat také tak, že si nejdříve vytvoříme komplexní proměnnou s: s = tf('s')

Přenos pak zadáme následujícím způsobem: sys=(1.2*s+0.5)/(7.1*s^3+2*s^2+1.5*s+0.1) Transfer function: 1.2 s + 0.5 ----------------------------7.1 s^3 + 2 s^2 + 1.5 s + 0.1

Zadání diskrétního přenosu za pomocí funkce tf Zadání diskrétního systému je obdobné jako u systému spojitého. Navíc zde zadáváme vzorkovací periodu Ts (častěji bývá vzorkovací perioda označována T v matlabu je však uvedeno značení Ts). Základní syntaxe příkazu: SYS = TF(NUM,DEN,Ts)

– diskrétní systém, zadáváme čitatele jmenovatele a

vzorkovací periodu z = TF('z',Ts) – vytvoření diskrétní proměnné z

Analýza systémů

20

Máme systém zadaný přenosem: G( z) =

1,2 z + 0,5 7,1z + 2 z 2 + 1,5 z + 0,1 3

Pro zadání přenosu použijeme příkaz: sys = tf([1.2 0.5], [7.1 2 1.5 0.1], 0.1)

Do proměnné sys se uložil přenos v tomto tvaru: Transfer function: 1.2 z + 0.5 ----------------------------7.1 z^3 + 2 z^2 + 1.5 z + 0.1 Sampling time: 0.1

Zadání přenosu v programu Simulink Funkce tf je v programu Simulink zprostředkována pomocí bloků Transfer Fnc pro spojité systémy a Disckrete Transfer Fnc pro diskrétní systémy viz. obr. 3.11.

obr. 3.11 Bloky tf programu Simulink

Na obrázku obr. 3.12 je příklad zapojení bloku Transfer Fnc pro vykreslení přechodové charakteristiky spojitého systému.

obr. 3.12 Blokové schéma pro vykresleni přechodové charakteristiky pomocí bloku Transfer Fnc

3.1.3 Zadání přenosu systému pomocí pólů a nul Je-li přenos systému vyjádřen tvarem G (s) =

bm ( s + n1 ) ... ( s + nm ) , an ( s + p1 ) ... ( s + pn )

(3.14)

což je rozklad čitatele i jmenovatele na kořenové činitele, vidíme z něj rozložení pólů a nul. Toto rozložení zcela charakterizuje dynamické vlastnosti systému. Nuly a póly mohou být reálné, komplexní nebo i ryze imaginární. Pokud jsou komplexní a ryze imaginární, musí být vždy po dvou komplexně sdružené. V komplexní rovině označujeme nuly kroužkem a póly křížkem. Čím jsou póly dále od imaginární osy, tím je přechodový děj více tlumen a tím kratší. Jsou-li póly komplexní znamená to, že přechodový děj bude mít kmitavou složku. Jsou-li nuly blíže imaginární osy než póly, bude převládat derivační složka. Póly v pravé polovině znamenají vždy nestabilní pochod. Póly v počátku vyjadřují integrační charakter, nuly v počátku pak zase derivační charakter. [Švarc, I. 1992]

Analýza systémů

21

Příklady: Zadání spojitého přenosu za pomocí funkce zpk - create zero-pole-gain model. SYS = ZPK(Z,P,K) – spojitý systém, zadáváme nuly, s = ZPK('s') – vytvoření Laplaceovy proměnné s

póly a zesílení

Máme systém zadaný přenosem: G (s) =

3s 2 ( s + 2) s 2 ( s + 5)( s + 6)( s + 2)

Pro zadání přenosu použijeme příkaz: sys=zpk([0 0 -2],[0 0 -5 -6 -2],3) Zero/pole/gain: 3 s^2 (s+2) --------------------s^2 (s+5) (s+6) (s+2)

Při zadávání přenosu pomocí nul, pólů a zesílení můžeme nejdřív zadat do vlastních proměnných hodnoty nul, pólů a zesílení pomocí těrcho proměnných potom zadáme přenos: z=[0 0 -3]; p=[0 0 -5 -6 -2]; k = 3; sys=zpk(z,p,k) Zero/pole/gain: 3 s^2 (s+3) --------------------s^2 (s+5) (s+6) (s+2)

Přenos můžeme zadat také tak, že si nejdříve vytvoříme Laplaceovu proměnnou s: s = zpk('s') Zero/pole/gain: s

Přenos pak zadáme následujícím způsobem: H = (3*s^2*(s+3))/(s^2*(s+5)*(s+6)*(s+2)) Zero/pole/gain: 3 s^2 (s+3) --------------------s^2 (s+5) (s+6) (s+2)

Zadání diskrétního přenosu za pomocí funkce zpk - create zero-pole-gain model. SYS = ZPK(Z,P,K,Ts)

– diskrétní systém, zadáváme nuly, póly, zesílení

a vzorkovací periodu z = ZPK('z',Ts) – vytvoření diskrétní Laplaceovy proměnné z Máme systém zadaný přenosem: G( z) =

3 z 2 ( z + 2) z 2 ( z + 5)( z + 6)( z + 2)

Pro zadání přenosu použijeme příkaz: sys=zpk([0 0 -2],[0 0 -5 -6 -2],3,0.1)

Analýza systémů

22

Výsledkem je přenos ve tvaru: Zero/pole/gain: 3 z^2 (z+2) --------------------z^2 (z+5) (z+6) (z+2) Sampling time: 0.1

Zadání přenosu pomocí funkce zpk v programu Simulink Funkce zpk je v programu Simulink zprostředkována pomocí bloků Zero-Pole pro spojité systémy a Diskrete Zero-Pole pro systémy diskrétní viz. obr. 3.13.

obr. 3.13 Bloky Zero-Pole programu Simulink

Na obrázku obr. 3.14 je příklad zapojení bloku Zero-Pole pro vykreslení přechodové charakteristiky spojitého systému.

obr. 3.14 Blokové schéma zapojení bloku Zero-Pole pro vykreslení přechodové charakteristiky

3.1.4 Zadání stavového popisu Stavový popis spojitého systému je plně určen čtveřicí matic A,B,C,D. V běžných případech je matice D=0. Pro zápis stavového popisu používáme soustavu rovnice ve tvaru (3.15). Matice A je matice systému, matice B je matice buzení, C je matice výstupu a D matice převodu. dx = Ax + Bu dt y = Cx + Du .

x& =

(3.15)

Příklady: Zadání stavového popisu pro spojitý systém za pomocí funkce ss - Create statespace model. SYS = SS(A,B,C,D) -

zadání spojitého systému pomocí matic A,B,C,D

Máme systém zadaný následující soustavou rovnic: x&1 (t ) = −2 x1 (t ) − x2 (t ) + u (t ) x& 2 (t ) = x1 (t ) y (t ) = x1 (t ) + 2 x2 (t )

Analýza systémů

23

Nejdříve do proměnných A,B,C,D Zadáme hodnoty jednotlivých matic: A= [-2, -1; 1, 0]; B = [1; 0]; C = [1, 2]; D = [0];

Nakonec pomocí těchto matic zadáme systém: SYS = SS (A, B, C, D) a = x1 x2

x1 -2 1

x1 x2

u1 1 0

y1

x1 1

x2 -1 0

b =

c = x2 2

d = u1 y1 0 Continuous-time model.

Zadání stavového popisu pro diskrétní systém za pomocí funkce ss - Create statespace model. SYS = SS(A,B,C,D,Ts) -

zadání diskrétního systému pomocí matic A,B,C,D a

vzorkovací periody Ts. Máme systém zadaný následující soustavou rovnic: x&1 (k + 1) = −2 x1 (k ) − x2 (k ) + u (k ) x& 2 (k + 1) = x1 (k ) y (k ) = x1 (k ) + 2 x2 (k ) Nejdříve do proměnných A,B,C,D Zadáme hodnoty jednotlivých matic: A= [-2, -1; 1, 0]; B = [1; 0]; C = [1, 2]; D = [0];

Analýza systémů

24

Nakonec pomocí těchto matic a hodnoty Ts zadáme systém: SYS = SS (A, B, C, D, 0.1) a = x1 x2

x1 -2 1

x1 x2

u1 1 0

y1

x1 1

x2 -1 0

b =

c = x2 2

d = u1 y1 0 Sampling time: 0.1 Discrete-time model.

Zadání stavového popisu systému v programu Simulink Funkce ss je v programu Simulink zprostředkována pomocí bloků State-Space a Discrete State-Space viz.obr. 3.15.

obr. 3.15 Bloky State-Space programu Simulink

Na obrázku obr. 3.16 je příklad zapojení bloku State-Space pro vykreslení přechodové charakteristiky spojitého systému.

obr. 3.16 Blokové schéma zapojení bloku State-Space pro vykreslení přechodové charakteristiky

3.1.5 Převody mezi jednotlivými způsoby popisu systému Máme-li již zadaný popis systému pomocí jednoho z výše uvedených způsobu (tf, zpk, ss), a potřebujeme-li jej vyjádřit i v jiném tvaru, můžeme využít funkcí programu Matlab, které tento převod umožňují. Dále můžeme převádět spojité systémy na diskrétní a naopak. Diskrétní systémy pak můžeme převádět na systémy s jinou vzorkovací periodou.

Převod mezi stavovým popisem a přenosem zadaným pomocí funkce tf Nejdříve zadáme stavový popis který budeme chtít převést na přenos: A=[-1 -3 -2;1 0 0;0 1 0]; B=[1; 0; 0]; C=[0 4 7]; D=0; SSsys = ss(A,B,C,D)

Jednou z možností jak stavový popis převést je použití funkce ss2tf. Pokud použijeme tuto funkci musíme nejdřív znát všechny matice stavového popisu. Pro jejich získání můžeme použít příkaz ssdata: [A,B,C,D]=ssdata(SSsys)

Analýza systémů

25

V dalším kroku pak použijeme příkaz ss2tf pro získání čitatele a jmenovatele přenosu [NUM, DEN] = ss2tf (A, B, C, D)

posledním krokem je zadání přenosu pomocí funkce tf TFsys=tf(NUM, DEN) Transfer function: 4 s + 7 ------------------s^3 + s^2 + 3 s + 2

Druhou možností je převést stavoví popis na přenos přímo pomocí funkce tf: TFsys = tf(SSsys)

Převod mezi přenosem zadaným pomocí funkce tf a stavovým popisem Začneme zadáním přenosu systému: TFsys = tf([2 5], [7 2 5 1])

Při převodu pomocí funkce tf2ss je postup následující: Pomocí funkce tfdata si vyjádříme čitatele a jmenovatele přenosu (parametr v ukládá hodnoty jako vektory) [NUM, DEN]=tfdata(TFsys, 'v')

pak pomocí funkce tf2ss získáme matice stavového popisu A, B, C, D [A,B,C,D] = tf2ss (NUM,DEN)

nakonec zadáme stavový popis pomocí příkazu ss: SSsys=ss(A,B,C,D)

Druhou možností je převod přímo pomocí funkce ss: SSsys = ss(TFsys) a = x1 x2 x3

x1 -0.28571 1 0

x1 x2 x3

u1 1 0 0

y1

x1 0

x2 -0.71429 0 1

x3 -0.14286 0 0

x2 0.28571

x3 0.71429

b =

c =

d = y1 Continuous-time model.

u1 0

Převod mezi stavovým popisem a přenosem zadaným pomocí funkce zpk Postup je obdobný jako v předcházejících případech: Zadání stavového popisu A=[-1 -3 -2;1 0 0;0 1 0]; B=[1; 0; 0]; C=[0 4 7]; D=0; SSsys = ss(A,B,C,D)

Analýza systémů

26

Převod pomocí funkce ss2zp. V druhém kroku se do proměnných Z,P,K ukládají hodnoty nul, pólů a zesílení. [A,B,C,D]=ssdata(SSsys) [Z,P,K] = ss2zp (A, B, C, D) ZPKsys=zpk(Z,P,K)

Převod pomocí funkce zpk ZPKsys1 = zpk(SSsys)

Převod mezi přenosem zadaným pomocí funkce zpk a stavovým popisem Zadání přenosu pomocí funkce zpk ZPKsys = zpk([0 0 -2],[0 0 -5 -6 -2],3)

Převod pomocí funkce zp2ss [Z,P,K]=zpkdata(ZPKsys,'v') [A,B,C,D]=zp2ss(Z,P,K) SSsys=ss(A,B,C,D)

Převod pomocí funkce ss SSsys1 = ss(ZPKsys)

Převod mezi přenosem zadaným pomocí funkce tf a přenosem z funkce zpk Zadání přenosu pomocí funkce tf TFsys = tf([2 5], [7 2 5 1])

Převod pomocí funkce tf2zp [NUM, DEN]=tfdata(TFsys,'v') [Z,P,K]=tf2zp(NUM, DEN) ZPKsys=zpk(Z,P,K)

Převod pomocí funkce zpk ZPKsys1 = zpk(TFsys)

Převod mezi přenosem zadaným pomocí funkce zpk a přenosem z funkce tf Zadání přenosu pomocí funkce zpk ZPKsys=zpk([0 0 -2],[0 0 -5 -6 -2],3)

Převod pomocí funkce zp2tf [Z,P,K]=zpkdata(ZPKsys,'v') [NUM, DEN]=zp2tf(Z,P,K) TFsys=tf(NUM, DEN)

Převod pomocí funkce tf TFsys1 = tf(ZPKsys)

Převod spojitého systému na diskrétní Pro převod spojitého systému na diskrétní je určena funkce c2d. Příkaz zadáváme ve tvaru SYSD = c2d(SYSC,Ts),

kde SYSD je výsledný diskrétní přenos, SYSC je zdrojový přenos, Ts vzorkovací frekvence.

Analýza systémů

27

Příklad: Převod spojitého systému na diskrétní Máme spojitý systém zadaný přenosem ve tvaru: G (s) =

2s + 5 . 7 s + 2 s 2 + 5s + 1 3

Pro zadání přenosu použijeme funkci tf: sysC = tf([2 5], [7 2 5 1]) Transfer function: 2 s + 5 ----------------------7 s^3 + 2 s^2 + 5 s + 1

Převod na diskrétní tvar provedeme pomocí příkazu c2d se zadanou vzorkovací periodou 0,5: sysD = c2d(sysC,0.5) Transfer function: 0.0478 z^2 + 0.05293 z - 0.01876 ---------------------------------z^3 - 2.695 z^2 + 2.578 z - 0.8669 Sampling time: 0.5

Porovnání přechodové charakteristiky původního spojitého systému a získaného diskrétního systému (obr. 3.1) můžeme vykreslit pomocí příkazu step: step(sysC, sysD) Step Response From: U(1) 6

5

To: Y(1)

Amplitude

4

3

2

1

0 0

10

20

30

40

50

60

Time (sec.)

obr. 3.17 Porovnání přechodové charakteristiky spojitého a diskrétního přenosu

Převod diskrétního systému na spojitý Pro převod diskrétního systému na spojitý je určena funkce d2c. Příkaz zadáváme ve tvaru SYSC = d2c(SYSD),7

kde SYSC je výsledný diskrétní přenos a SYSD je zdrojový přenos.

Analýza systémů

28

Změna vzorkovací frekvence u diskrétního systému Pro změnu vzorkovací frekvence u diskrétních systémů je určen příkaz d2d. Příkaz zadáváme ve tvaru SYS = d2d(SYS, Ts),

kde SYS je upravovaný systém a Ts je nová vzorkovací frekvence.

3.2 Propojení systémů Složitější systémy není možné jednoduše zadat jedním přenosem či stavovým modelem. U složitějších úloh je třeba propojit několik systémů dohromady. Program Matlab nabízí několik možností jak jednotlivé systémy vzájemně spojit.

3.2.1 Paralelní zapojení Pro paralelní zapojení (obr. 3.19) jednotlivých systémů slouží v programu Matlab funkce parallel. U jednorozměrných modelů je možné místo funkce parallel použít součet. Paralelní zapojení dvou jednorozměrných systémů SYS1 a SYS2 provedeme pomocí příkazů: SYS = PARALLEL(SYS1,SYS2),

nebo SYS = SYS2 + SYS1

výsledný systém se uloží do proměnné SYS. SYS 1 u

+ y

SYS 2

+

obr. 3.18 Paralelní zapojení jednorozměrných systémů

Matlab pomocí nadstavby Control System umožňuje definovat a propojovat mnohorozměrné systémy. Příklad takového paralelního zapojení je vidět na obr. 3.19. v1 z1 SYS 1 y1 u1 + u y u2 + y2 v2 SYS 2 z2 obr. 3.19 Paralelní zapojení mnohorozměrných systémů

Analýza systémů

29

3.2.2 Sériové zapojení Pro sérové zapojení (obr. 3.20) jednotlivých systémů slouží v programu Matlab funkce series. U jednorozměrných modelů je možné místo funkce series použít násobení. Sériové zapojení dvou jednorozměrných systémů SYS1 a SYS2 provedeme pomocí příkazů: SYS = SERIES(SYS1,SYS2),

nebo SYS = SYS2 * SYS1

výsledný systém se uloží do proměnné SYS.

u

SYS 1

SYS 2

y

obr. 3.20 Sériové zapojení jednorozměrných systémů

Matlab pomocí nadstavby Control System umožňuje definovat a propojovat mnohorozměrné systémy. Příklad takového sériového zapojení je vidět na obr. 3.21. v2 y1 u2 SYS 2 y u SYS 1 z1 obr. 3.21 Sériové zapojení mnohorozměrných systémů

3.2.3 Zpětnovazební zapojení Pro zpětnovazební zapojení jednotlivých systémů slouží v programu Matlab funkce feedback. Na rozdíl od paralelního a sériového zapojení se nedá nahradit matematickým operátorem. Můžeme použít dva základní typy zpětnovazebního zapojení. Zapojení se zápornou vazbou (obr. 3.22) a zapojení s kladnou vazbou (obr. 3.23). Zpětnovazební zapojení se zápornou vazbou pro jeden systém SYS1 (obr. 3.22) deklarujeme následujícím příkazem: SYS = FEEDBACK(SYS1,1,-1).

Parametr –1 není třeba zadávat. Záporná vazba je nejpoužívanější, a proto je tato hodnota nastavena automaticky jako výchozí, aniž bychom ji zadali. + y u SYS 1 y _ obr. 3.22 Zpětnovazební zapojení systému se zápornou vazbou

Analýza systémů

30

Zapojení s kladnou zpětnou vazbou pro jeden systém SYS1 (obr. 3.23) deklarujeme následujícím příkazem: SYS = FEEDBACK(SYS1,1,+1)

u

+

SYS 1 +

y y

obr. 3.23 Zpětnovazební zapojení systému s kladnou vazbou Zpětnovazební zapojení se zápornou vazbou pro jeden systém SYS1 umístěný ve zpětné vazbě (obr. 3.24) deklarujeme tímto příkazem: SYS = FEEDBACK(1,SYS1)

u

+

y y

_ SYS 1

obr. 3.24 Zpětnovazební zapojení se systémem ve zpětné vazbě

Zpětnovazební zapojení se zápornou vazbou pro dva systémy SYS1 a SYS2 z nich jeden je umístěn ve zpětné vazbě (obr. 3.25) deklarujeme následujícím příkazem: SYS = FEEDBACK(SYS1,SYS2)

u

+

SYS 1

y

_ SYS 2 obr. 3.25 Zpětnovazební zapojení jednorozměrných systémů

Matlab pomocí nadstavby Control System umožňuje definovat a propojovat mnohorozměrné systémy. Příklad takového zpětnovazebního zapojení je vidět na obr. 3.26. v z SYS 1 + y u y _ SYS 2 obr. 3.26 Zpětnovazební zapojení mnohorozměrného systému

Analýza systémů

31

3.3 Vykreslení základních charakteristik systému V následující kapitole budou popsány možnosti vykreslení základních charakteristik systému: •

Přechodové charakteristiky

•

Impulsní charakteristiky

•

Kmitočtové charakteristiky

3.3.1 Přechodová charakteristika Přechodová charakteristika je odezva systému na jednotkový Heavisideův skok na vstupu. Je to tedy časový průběh výstupní veličiny systému a značí se h(t). Jednotkový skok je funkce, která do času t=0 má nulovou hodnotu a v čase t=0 skočí její hodnota na jednotku, kterou pak stále udržuje. [Švarc, I. 1992] V programu Matlab můžeme přechodovou charakteristiku vykreslit příkazem step, nebo pomocí příkazu ltiview s parametrem step. Tyto příkazy zadáváme v následujících tvarech: step(SYS) ltiview(‘step’,SYS)

Příklady: Proporcionální systém Provedeme porovnání přechodových charakteristik u tří proporcionálních systémů, které máme zadány přenosy: 1. proporcionální člen bez setrvačnosti (ideální) Gp( s ) = 0,9 2. proporcionální člen se setrvačností 1. řádu Gp1 ( s ) =

1 3s + 1

3. proporcionální člen se setrvačností 2. řádu Gp2 ( s ) =

5s + 1 7 s + 3s + 5 2

Pro zadání všech tří přenosů použijme příkaz tf : p = tf(0.9, 1) p1 = tf([1], [3 1]) p2 = tf([5 1], [7 3 5])

Přechodové charakteristiky vykreslíme pomocí příkazu step. step(p, p1, p2)

Pomocí příkazů title, xlabel a ylabel můžeme do grafu zadat titulek a popis osy x a y. title('Prechodova charakteristika') ylabel('h(t)') xlabel('cas[s]')

Analýza systémů

32

Výsledný graf je na obr. 3.27, popisky jednotlivých charakteristik byly přidány dodatečně pomocí příkazů z panelu nástrojů dialogového okna grafu. Pomocí dalších nabídek tohoto dialogového okna je možné měnit i některé další zobrazené informace např. zobrazení mřížky, popisky os a titulku, barvy vykreslených grafů a další. Prechodova charakteristika From: U(1) 1.2

1

Gp(s)

0.8

h(t)

To: Y(1)

Gp1(s) 0.6

Gp2(s) 0.4

0.2

0

-0.2 0

5

10

15

20

25

30

cas[s]

obr. 3.27 Vykreslení přechodových charakteristik proporcionálních systémů

Pro vykreslení přechodových charakteristik pomocí programu Simulink budeme potřebovat bloky Transfer Fnc pro každý přenos, generátor vstupního signálu Step, pro spojení výsledných dat blok Mux a pro zobrazení dat blok Scope. Do bloků Transfer Fnc zadáme hodnoty jednotlivých přenosů a bloky propojíme podle schématu na obr. 3.28.

obr. 3.28 Blokově schéma programu Simulink - vykresleni h(t) pro proporcionální člen

Integrační systém Provedeme porovnání přechodových charakteristik u tří integračních systémů, jež máme zadány přenosy: 1. Integrační člen bez setrvačnosti (ideální): Gi ( s ) = 2. Integrační člen se setrvačností 1. řádu: Gi1 ( s ) =

1 s

1 s (s + 1)

3. Integrační člen 2. řádu se setrvačností 1.řádu: Gi2 ( s ) =

1 s (s + 1) 2

Přenosy zadáme pomocí příkazu tf: i = tf([1],[1 0]); i1 = tf([1], [1 1 0]); i2 = tf([1], [1 1 0 0])

Analýza systémů

33

Pro vykreslení charakteristik použijeme příkaz step: step(i, i1, i2)

Výsledný graf ke kterému byly přidány popisky os a titulek je na obr. 3.29. Prechodova charakteristika From: U(1) 14

12

To: Y(1)

h(t)

10

8

Gi2(s)

6

Gi(s)

4

Gi1(s)

2

0 0

1

2

3

4

5

cas[s]

obr. 3.29 Vykreslení přechodových charakteristik integračních systémů

Pro vykreslení přechodových charakteristik pomocí programu Simulink budeme potřebovat bloky Transfer Fnc, Step, Mux a Scope. Bloky propojíme podle schématu na obr. 3.30.

obr. 3.30 Blokově schéma programu Simulink - vykresleni h(t) pro integrační člen

Derivační systém Podobně jako v předchozích případech porovnáme charakteristiky tří systémů tentokrát derivačních, které jsou zadány přenosy: 1. Derivační člen 1. řádu se setrvačností 1. řádu: Gd ( s ) =

3s 6s + 5

2. Derivační člen 1.řádu se setrvačností 2.řádu: Gd1 ( s ) =

3s 6 s + 3s + 1

3. Derivační člen se setrvačností 2. řádu: Gd 2 ( s ) =

2

3s + 1s 6 s + 3s 2 + 5s + 1 3

Stejně jako v předchozích případech použijeme pro zadání přenosu příkaz tf: d=tf([3 0],[6 3]); d1=tf([3 0],[6 3 1]);d2=tf([3 1 0],[6 3 5 1])

a pro vykreslení grafů příkaz step: step(d, d1, d2)

Výsledný graf je na obr. 3.31.

Analýza systémů

34 Prechodova charakteristika From: U(1) 0.7 0.6

Gd1(s) 0.5

Gd2(s)

0.4

Gd(s)

h(t)

To: Y(1)

0.3 0.2 0.1 0 -0.1 -0.2 -0.3 0

5

10

15

20

25

30

35

40

cas[s]

obr. 3.31 Vykreslení přechodových charakteristik derivačních systémů

V programu Simulink provedeme zadání podle schématu na obr. 3.32.

obr. 3.32 Blokové schéma programu Simulink - vykresleni h(t) pro deriva ční člen

3.3.2 Impulsní charakteristika Impulsní charakteristika g(t) je odezva systému na Diracův impuls na vstupu. Diracův nebo též jednotkový impuls se nedá fyzikálně realizovat. Při jeho realizaci se udává, že musí mít co největší amplitudu a nejkratší dobu trvání. Matematicky dosahuje Diracův impuls nenulové hodnoty pouze v čase t=0. [Švarc, I. 1992] V programu Matlab můžeme impulsní charakteristiku zobrazit pomocí příkazu impulse, nebo pomocí příkazu ltiview s parametrem impulse. Vykreslení provedeme pomocí příkazů: impulse(SYS) ltiview(‘impulse’,SYS)

Příklady: Proporcionální systém Podobně jako při vykreslení přechodové charakteristiky porovnáme charakteristiky tří proporcionálních systémů. Použijeme stejné přenosy, které zadáme pomocí funkce tf: Gp( s ) = 0,9 ; Gp1 ( s ) =

1 5s + 1 ; Gp2 ( s ) = 2 , 3s + 1 7 s + 3s + 5

Analýza systémů

35

pro jejich zadání použijeme opět funkci tf: p = tf(0.9, 1); p1 = tf([1], [3 1]); p2 = tf([5 1], [7 3 5])

Vykreslení pak provedeme pomocí příkazu impulse: impulse(p,p1,p2)

Výsledek podobě v grafu je na obr. 3.33. Impulsni charakteristika From: U(1) 1

0.8

0.6

g(t)

To: Y(1)

Gp2(s) 0.4

Gp1(s) 0.2

0

Gp(s) -0.2

-0.4 0

5

10

15

20

25

30

cas[s]

obr. 3.33 Impulsní charakteristiky proporcionálních systémů

Integrační systém Porovnání zadaných přenosů: Gi( s ) =

1 1 1 ; Gi2 ( s ) = 2 ; Gi1 ( s ) = s s (s + 1) s (s + 1)

Přenosy zadáme pomocí příkazu tf: i = tf([1],[1 0]); i1 = tf([1], [1 1 0]); i2 = tf([1], [1 1 0 0])

Vykreslení pak provedeme pomocí příkazu impulse: impulse(i,i1,i2)

Výsledek podobě v grafu je na obr. 3.34. Impulsni charakteristika From: U(1) 4.5 4 3.5

To: Y(1)

g(t)

3 2.5 2

Gi2(s)

Gi(s)

1.5 1 0.5

Gi1(s)

0 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

cas[s]

obr. 3.34 Impulsní charakteristiky integračních systémů

Analýza systémů

36

Derivační systém Gd ( s ) =

3s 3s 3s + 1s ; Gd1 ( s ) = 2 ; Gd 2 ( s ) = 3 6s + 5 6 s + 3s + 1 6 s + 3s 2 + 5s + 1

Stejně jako v předchozích případech použijeme pro zadání přenosu příkaz tf: d=tf([3 0],[6 3]); d1=tf([3 0],[6 3 1]);d2=tf([3 1 0],[6 3 5 1])

Grafy vykreslíme funkcí impulse: impulse(d, d1,d2)

Impulsni charakteristika From: U(1) 0.5 0.4 0.3

Gd1(s) Gd2(s)

g(t)

To: Y(1)

0.2 0.1 0 -0.1 -0.2 -0.3

Gd(s)

-0.4 0

5

10

15

20

25

30

35

40

cas[s]

obr. 3.35 Impulsní charakteristiky derivačních systémů

3.3.3 Kmitočtové charakteristiky Kmitočtový přenos je poměr výstupních a vstupních harmonických kmitů systému. G ( jω ) =

Y ( jω ) U ( jω )

Základním vstupním signálem, kterého používáme při kmitočtovém vyšetřování vlastností lineárních systémů je sinusový průběh. Vzhledem k lineárnosti uvažovaného systému je výstupní signál také sinusový a mění se jen jeho amplitudová fáze. Je-li vstupní signál u = u0 sin ωt je výstupní signál y = y0 sin(ωt + ϕ ) , kde u0, y0 jsou amplitudy vstupního a výstupního signálu, ω úhlová frekvence a ϕ fázové posunutí výstupního signálu vůči vstupnímu.

Analýza systémů

37

Kmitočtová charakteristika Kmitočtová charakteristika (obr. 3.36) je grafické znázornění G(jω) v komplexní (Gausově) rovině, přičemž proměnným parametrem je úhlová frekvence ω, kterou měním od ω=0 do ω=∞. Při sestrojování kmitočtové charakteristiky si kmitočtový přenos G(jω) upravíme na složkový tvar komplexního čísla G ( jω ) = Re(ω ) + j Im(ω ) , kde vhodným hodnotám ω počítáme reálnou a které zakreslujeme do komplexní roviny. [Švarc, I. 1992]

imaginární

složku

přenosu,

Re ω=∞

ω=0

Im

obr. 3.36 Kmitočtová charakteristika v komplexní rovině

Pro vykreslení kmitočtové charakteristiky v programu Matlab můžeme použit příkaz nyquist, nebo příkaz ltiview s parametrem nyquist. nyquist(SYS) ltiview(‘nyquist’,SYS)

Logaritmické kmitočtové charakteristiky Kmitočtovou charakteristiku také někdy nazýváme amplitudo-fázovou charakteristikou v komplexní rovině. Teoreticky i prakticky se dá rozdělit na dvě charakteristiky – amplitudovou a fázovou. Pro větší přehlednost tuto dvojici charakteristiky vykreslujeme v logaritmických souřadnicích (obr. 3.38). U obou charakteristik vynášíme hodnotu ω na vodorovnou osu v logaritmickém měřítku jako logω. Tím dosáhneme velké rozmezí změn frekvence. Na svislou osu amplitudové logaritmické charakteristiky vynášíme absolutní hodnotu kmitočtového přenosu, ale v jednotkách decibel. Decibel je v tomto případě definován jako dvacetinásobek dekadického logaritmu. A(ω ) = 20 log G ( jω ) Hodnoty A[dB] vynášíme již v lineárním měřítku. Fázi vynášíme na svislou osu opět v lineárním měřítku. [Švarc, I. 1992] Pro vykreslení logaritmické kmitočtové charakteristiky v programu Matlab můžeme použit příkaz bode, nebo příkaz ltiview s parametrem bode: bode(SYS) ltiview(‘bode’,SYS)

Analýza systémů

38

Příklady: Proporcionální systém Opět použijeme přenosy z předchozích kapitol, které zadáme pomocí funkce tf: Gp( s ) = 0,9 ; Gp1 ( s ) =

1 5s + 1 ; Gp2 ( s ) = 2 3s + 1 7 s + 3s + 5

p = tf(0.9, 1); p1 = tf([1], [3 1]); p2 = tf([5 1], [7 3 5])

Vykreslení kmitočtové charaktericky (obr. 3.37) pak provedeme pomocí příkazu nyquist: nyquist(p, p1, p2)

Logaritmicko kmitočtovou charakteristiku (obr. 3.38) vykreslíme pomocí funkce bode: bode(p,p1,p2) Kmitoctova charakteristika From: U(1) 1.5

1

Im

To: Y(1)

0.5

0

-0.5

-1

-1.5 -1

-0.5

0

0.5

1

1.5

2

Re

obr. 3.37 Kmitočtová charakteristika proporcionálních systémů Logaritmicko kmitoctova charakteristika 10 0

Faze (stupne); A (dB)

-10 -20 -30 50

0

-50

-100 10-2

10-1

100

101

Frekvence (rad/sec)

obr. 3.38 Logaritmická kmitočtová charakteristika proporcionálních systémů

Analýza systémů

39

Integrační systém Gi( s ) =

1 1 1 ; Gi1 ( s ) = ; Gi2 ( s ) = 2 s s (s + 1) s (s + 1)

Přenosy zadáme pomocí příkazu tf: i = tf([1],[1 0]); i1 = tf([1], [1 1 0]); i2 = tf([1], [1 1 0 0])

Vykreslení kmitočtové charaktericky pak provedeme pomocí příkazu nyquist: nyquist(i, i1, i2)

Logaritmickou kmitočtové charakteristiku vykreslíme pomocí funkce bode: bode(i, i1, i2)

Derivační systém Pro porovnání máme opět zadány následující systémy: Gd ( s ) =

3s 3s 3s + 1s ; Gd1 ( s ) = 2 ; Gd 2 ( s ) = 3 6s + 5 6 s + 3s + 1 6 s + 3s 2 + 5s + 1

Stejně jako v předchozích případech použijeme pro zadání přenosu příkaz tf: d=tf([3 0],[6 3]); d1=tf([3 0],[6 3 1]);d2=tf([3 1 0],[6 3 5 1])

Vykreslení kmitočtové charaktericky pak provedeme pomocí příkazu nyquist: nyquist(i, i1, i2)

Logaritmicko kmitočtovou charakteristiku vykreslíme pomocí funkce bode: bode(i, i1, i2)

3.4 Stabilita systému Stabilita je vlastnost systému. Je jedním ze základních požadavků, které klademe na regulační obvod. Regulační obvod je stabilní, jestliže po vychýlení regulačního obvodu z rovnovážného stavu a odstranění podnětu, který tuto odchylku způsobil, se regulační obvod během času vrátí do původního rovnovážného stavu. Jinak řečeno je stabilita vlastnost regulačního obvodu vrátit se do rovnovážného stavu, jestliže skončí působení vzruchu, který ho z rovnovážného stavu vyvedl. [Švarc, I. 1992] Podle stability rozlišujeme tři typy systémů: stabilní, na mezi stability a nestabilní. Nutná a postačující podmínka stability pro spojité systémy: regulační obvod je stabilní, leží-li všechny kořey charakteristického polynomu v záporné polovině komplexní roviny (obr. 3.39). Im

0

Re

obr. 3.39 Oblast stabilních kořenů pro spojité systémy

Analýza systémů

40

Nutná a postačující podmínka stability pro diskrétní systémy: systém je stabilní právě tehdy, když jsou kořeny charakteristického mnohočlenu v absolutní hodnotě menší než jedna (obr. 3.40). 1

Im

Re 0

-1

1

-1

obr. 3.40 Oblast stabilních kořenů pro diskrétní systémy

Pro určení stability v programovém prostředí Matlab můžeme užít příkaz roots, pomocí něhož určíme hodnotu kořenů charakteristického polynomu. Na základě nutné a postačující podmínky stability, pak můžeme určit, zda jde o stabilní systém či nikoliv. Pomocí této podmínky určím pouze, zda je systém stabilní, pokud však není, nelze určit zda se jedná o systém na mezi stability nebo o systém nestabilní. Pro spojité systémy je další možností Nyquistovo kriterium stability. Jedná se o kmitočtové kriterium, které rozhoduje o stabilitě regulačního obvodu na základě průběhu kmitočtové charakteristiky otevřené smyčky obvodu. Podle Nyquistova kriteria musí pro stabilní obvod platit: 1.

2.

3.

Im{Go( jω )} = 0 ⇒ ω krit Re{Go(ω krit )} < −1

ϕ (ω krit ) = −π ⇒ ω krit A(ω krit ) < 1 ϕ (ω krit ) = −π ⇒ ω krit L(ω krit ) < 0

3.4.1 Stabilní systém Máme regulovanou soustavu zadanou přenosem GS(s), regulátor je typu PI zadaný přenosem GR(s) v obvodu je zapojen měřící člen s přenosem Gmč(s). Regulační obvod je zapojen podle obr. 3.41. G R ( s) =

1 5 1 ; GS (s) = ; Gmč ( s ) = 5s s +1 2s + 1

+

-

GR(s)

GS(s) Gmč(s)

obr. 3.41 Zapojení regulačního obvodu

Analýza systémů

41

Ověření pomocí v programem Matlab Nejdříve pomočí funkce tf zadáme jednotlivé přenosy: R = tf([1],[5 0]); S = tf([5],[1 1]); mc =tf([1],[2 1]);

Jednotlivé přenosy propojíme podle zadání a získáme přenos regulačního obvodu. RS = series(R,S); sys = feedback (RS,mc)

Určení stability pomocí kořenů charakteristického polynomu Nejdříve určíme stabilitu pomocí nutné a postačující podmínky stability, tedy pomocí kořenů charakteristického polynomu. Pomocí příkazu tfdata určíme hodnotu charakteristického polynomu, která se nám načte do proměnné jmen. [cit,jmen]=tfdata(sys,'v');

Kořeny polynomu pak zjistíme pomocí příkazu roots: r = roots (jmen) r = -1.3982 -0.0509 + 0.5958i -0.0509 - 0.5958i

Z výsledku vidíme, že se jedná o stabilní systém, neboť všechny reálné složky kořenů mají zápornou hodnotu:

Určení stability pomocí Nyquistova kriteria Pro určení stability podle Nyquistova kriteria stability si nejdříve vypočítáme přenos otevřeného obvodu: Go=series(RS,mc)

Pomocí příkazu ltiview vykreslíme amplitudo fázovou kmitočtovou charakteristiku. ltiview('nyquist',Go)

Na výsledném grafu (obr. 3.42) vidíme že se jedná o stabilní systém neboť kmitočtová charakteristika protíná reálnou osu až za hodnotou –1.

Analýza systémů

42 Nyquist Diagrams From: U(1) 10

To: Y(1)

Imaginary Axis

5

0

-5

-10

-3

-2.5

-2

-1.5

-1

-0.5

0

Real Axis

obr. 3.42 Kmitočtová charakteristika stabilního systému

Stabilitu můžeme určit i z logaritmické kmitočtové charakteristiky vykreslené pomocí příkazu ltiview: ltiview('bode',Go)

Na výsledném grafu (obr. 3.43) je vidět, že hodnota zesílení pro ϕ=180° je menší než 0, jedná se tedy o stabilní systém. Bode Diagrams 50

Phase (deg); Magnitude (dB)

0

-50

-100 0

-100

-200

-300 10-2

10-1

100

101

Frequency (rad/sec)

obr. 3.43 Logaritmická kmitočtová charakteristika stabilního systému

Ověření v programu Simulink Pro ověření stability v programu Simulink použijme bloky Step, Transfer Fnc, Sum a Scope zapojené podle schéma na obr. 3.44.

Analýza systémů

43

obr. 3.44 Zapojení bloků v programu Simulink

Výsledkem je přechodová charakteristika, na níž je vidět, jak se systém po určité době ustálí. Jedná se tedy o stabilní systém. Prechodova charakteristika 2.5

2

h(t)

1.5

1

0.5

0 0

20

40

60 cas[s]

80

100

120

obr. 3.45 Přechodová charakteristika stabilního systému

3.4.2 Systém na mezi stability Pro určení stability v tomto příkladu použijeme stejné zadání jako v předchozím případě pouze změníme hodnotu regulátoru na TI = 10/3. G R ( s) =

1 5 1 ; GS (s) = ; Gmč ( s ) = 10 s +1 2s + 1 s 3

Ověření pomocí v programem Matlab Nejdříve zadáme přenos regulačního obvodu: %nastaveni regulatoru TI = 10/3; R = tf([1],[TI 0]); %zadání soustavy a mc S = tf([5],[1 1]); mc =tf([1],[2 1]); %propojeni RS = series(R,S); sys = feedback (RS,mc)

Analýza systémů

44

Určení stability pomocí kořenů charakteristického polynomu Použijeme stejný postup jako v předchozím příkladu: [cit,jmen]=tfdata(sys,'v'); r1 = roots (jmen) r1 = -1.5000 0.0000 + 0.7071i 0.0000 - 0.7071i

Z výsledku vidíme, že se jedná o nestabilní systém neboť pouze jeden kořen splňuje podmínku stability.

Určení stability pomocí Nyquistova kriteria Nejdříve si určíme přenos otevřené smyčky Go(s): Go=series(RS,mc);

Pomocí následujícího příkazu nejdříve vykreslíme kmitočtovou charakteristiku: ltiview('nyquist',Go)

Z kmitočtové charakteristiky (obr. 3.46) vidíme, že se jedná o systém na mezi stability neboť charakteristika protíná reálnou osu v bodě –1. Nyquist Diagrams From: U(1)

3

1 To: Y(1)

Imaginary Axis

2

0

-1

-2

-3 -3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

Real Axis

obr. 3.46 Kmitočtová charakteristika systému na mezi stability

Stabilitu lze určit i z logaritmické kmitočtové charakteristiky: ltiview('bode',Go)

Z výsledného grafu na obr. 3.47 vidíme, že jde o systém na mezi stability, neboť zesílení pro ϕ=180° je rovno nule.

Analýza systémů

45 Bode Diagrams From: U(1) 50

-50

-100 0

To: Y(1)

Phase (deg); Magnitude (dB)

0

-100

-180 -200

-300 10-2

10-1

100

101

Frequency (rad/sec)

obr. 3.47 Logaritmická kmitočtová charakteristika systému na mezi stability

Ověření v programu Simulink Zapojení bloků je stejné jako v předchozím příkladě pouze se změní hodnota parametrů regulátoru.

obr. 3.48 Zapojení bloků v programu Simulink

Na obr. 3.49 je vidět výsledná přechodová charakteristika systému na mezi stability, který kmitá s konstantní amplitudou.

Analýza systémů

46 Prechodova charakteristika

3 2.5 2

h(t)

1.5 1 0.5 0 -0.5 -1

0

10

20

30

40

50 cas[s]

60

70

80

90

100

obr. 3.49 Přechodová charakteristika systému na mezi stability

3.4.3 Nestabilní systém Opět použijeme podobné zadání jako v předchozích případech. Změníme však hodnotu parametrů regulátoru na TI=2. G R ( s) =

1 5 1 ; GS (s) = ; Gmč ( s ) = 2s s +1 2s + 1

Ověření pomocí v programem Matlab Nejdříve zadáme přenos jednotlivých prvků regulačního obvodu a propojíme je: TI = 2 R = tf([1],[TI 0]) S = tf([5],[1 1]); mc =tf([1],[2 1]); RS = series(R,S); %propojení regulátoru a soustavy sys = feedback (RS,mc) %zapojení zpětné vazby a měřícího členu

Analýza systémů

47

Určení stability pomocí kořenů charakteristického polynomu Použijeme stejný postup jako v předchozích příkladech: [cit,jmen]=tfdata(sys,'v'); r1 = roots (jmen) r1 = -1.6545 0.0772 + 0.8658i 0.0772 - 0.8658i

Z výsledku vidíme, že se jedná o nestabilní systém neboť dva kořeny mají reálné složky v kladné části komplexní roviny.

Určení stability pomocí Nyquistova kriteria Nejdříve si určíme přenos otevřené smyčky Go(s): Go=series(RS,mc);

Pomocí následujícího příkazu nejdříve vykreslíme kmitočtovou charakteristiku: ltiview('nyquist',Go)

Z výsledné kmitočtové charakteristiky (obr. 3.50) vidíme, že se jedná o nestabilní systém, protože charakteristika protíná reálnou osu před bodem –1. Nyquist Diagrams From: U(1) 8

6

2 To: Y(1)

Imaginary Axis

4

0

-2

-4

-6 -8 -7

-6

-5

-4

-3

-2

-1

0

Real Axis

obr. 3.50 Kmitočtová charakteristika nestabilního systému

Pokud budeme stabilitu určovat z logaritmické kmitočtové charakteristiky použijeme příkaz: ltiview('bode',Go)

Analýza systémů

48

Na výsledném grafu (obr. 3.47) vidíme, že jde o systém nestabilní, neboť logaritmická amplituda pro ϕ=180° je větší než nula. Bode Diagrams From: U(1) 50

-50

-100 0

To: Y(1)

Phase (deg); Magnitude (dB)

0

-100

180° -200

-300 10-2

10-1

100

101

Frequency (rad/sec)

obr. 3.51 Logaritmická kmitočtová charakteristika nestabilního systému

Ověření v programu Simulink Zapojení bloků na obr. 3.52 je stejné jako v předchozích příkladech pouze se změní hodnota parametrů regulátoru.

obr. 3.52 Zapojení bloků v programu Simulink

Na výsledné přechodové charakteristice (obr. 3.53) je dobře vidět, že se jedná o nestabilní systém, který se v čase neustálí, ale právě naopak se zvyšuje amplituda kmitů. Prechodova charakteristika 80 60 40 20

h(t)

0 -20 -40 -60 -80 -100 0

5

10

15

20

25 cas[s]

30

35

40

45

50

obr. 3.53 Přechodová charakteristika nestabilního systému

Syntéza regulačních obvodů

49

4 Syntéza regulačních obvodů 4.1 Návrh regulátoru metodou Ziegler – Nichols Jde o metodu založenou na analytickém nebo experimentálním určování kritického zesílení. Po zjištění kritického zesílení provedeme výpočet koeficientů spojitého i diskrétního regulátoru na základě tabulky 3.1. Metodu můžeme přímo v programu Matlab kombinovat s graficko-analytickým postupem, při kterém zjistíme nejdříve kritickou frekvenci a kritické zesílení na základě příkazů bode (logaritmická kmitočtová charakteristika) a margin (kritická hodnota zesílení a kritická hodnota frekvence), a potom z těchto hodnot určím koeficienty regulátoru podle tabulky 3.1 [Dušek F. 2000]. tab. 4.1 Nastavení koeficientů metodou Zieglera - Nicholse

Typ regulátoru

Proporcionální složka kp

Integrační složka TI

Derivační složka TD

P

0,50 kkrit

-

-

PI

0,45 kkrit

0,85 Tkrit

-

PD

0,40 kkrit

-

0,50 Tkrit

PID

0,60 kkrit

0,50 Tkrit

0,12 Tkrit

Kkrit – kritické zesílení; Tkrit – kritická perioda (Tkrit=2*π/ωkrit)

Příklad: Vypočítejte koeficienty spojitého PID regulátoru pro soustavu zadanou přenosem: Gs( s ) =

9 s ( s + 6 s + 9) 2

Řešení: Nejdříve zadáme přenos systému: sys=tf(9,[1 6 9 0]) Transfer function: 9 ----------------s^3 + 6 s^2 + 9 s

Pomocí příkazu margin určíme hodnoty kritického zesílení a kritické frekvence. [Kr,Fk,wk,wf]=margin(sys) Získali jsme Kr = 6 Fk = 56.0798 wk = 3 wf = 0.9149

následující hodnoty: % kritická hodnota zesílení % kritická fáze % kritické frekvence pro kritické zesílení % kritické frekvence pro kritickou fázi

Syntéza regulačních obvodů

50

Tyto hodnoty můžeme zobrazit i graficky pomocí příkazu: margin(sys) Bode Diagrams Gm=15.563 dB (at 3 rad/sec), Pm=56.08 deg. (at 0.91491 rad/sec) 50

Phase (deg); Magnitude (dB)

0 -50 -100 -150 -50 -100 -150 -200 -250 -300 10 -2

10 -1

100

101

10 2

Frequency (rad/sec)

obr. 4.1 Logaritmická kmitočtová charakteristika soustavy s vyznačením kritických hodnot

V dalším kroku vypočítáme kritickou periodu Tk. Tk=2*pi/wk

Nyní již známe obě potřebné hodnoty (kritické zesílení a perioda) pro výpočet parametrů regulátoru: K=0.6*Kr % vypocet zesileni PID regulatoru Ti=0.5*Tk % vypocet integracni casove konstanty Td=0.12*Tk % vypocet derivacni casove konstanty

Z vypočítaných hodnot k, TI, TD sestavíme přenos regulátoru: reg = tf([K*Ti*Td K*Ti K],[Ti 0]) % prenos PID regulatoru Transfer function: 0.9475 s^2 + 3.77 s + 3.6 ------------------------1.047 s

Posledním krokem je propojení soustavy a regulátoru do uzavřeného regulačního obvodu k tomu můžeme použít příkaz feedback pro zpětnovazební propojení. Go=sys*reg %prenos otevreneho regulacniho obvodu Gro=feedback(Go,1) %prenos uzavreneho regulacniho obvodu Transfer function: 8.527 s^2 + 33.93 s + 32.4 -------------------------------------------------1.047 s^4 + 6.283 s^3 + 17.95 s^2 + 33.93 s + 32.4

Přechodová charakteristika uzavřeného regulačního obvodu je na obr. 4.2.

Syntéza regulačních obvodů

51 Step Response From: U(1)

1.8 1.6 1.4

To: Y(1)

Amplitude

1.2 1 0.8 0.6 0.4 0.2 0 0

1.5

3

4.5

6

7.5

9

Time (sec.)

obr. 4.2 Přechodová charakteristika uzavřeného regulačního obvodu

Syntéza regulačních obvodů

52

4.2 Návrh regulátoru metodou požadovaného modelu Metoda umožňuje snadné a rychlé seřízení standardních typů číslicových a analogových regulátorů pro základní druhy regulovaných soustav s dopravním zpožděním. Typ regulátoru jde doporučen z hlediska vlastností regulované soustavy a požadavků na nulovou trvalou regulační odchylku způsobenou skokovou změnou žádané hodnoty. Předpokládá se použití standardních s odpovídajícími L- a Z- přenosy.

analogových

číslicových

a

regulátorů

Aby bylo možné metodu požadovaného modelu použít, musí být přenos regulované soustavy GS(s) v jednom ze základních tvarů uvedených v tab. 4.2. Pro některé přenosy které neodpovídají základním tvarům je možné určit přenos náhradní viz [Vítečková 1998]. tab. 4.2 Nastavení koeficientů regulátoru metodou inverzní dynamiky

Regulátor Regulovaná soustava

kp*

Typ Td=0

Td>0

P

2 k1 (2Tw + T )

a k1

k1 e −Td s (T1 s +1)

PI

2TI k1 (2Tw + T )

aTI k1

k1 e −Td s s (T1 s + 1)

PD

2 k1 (2Tw + T )

a k1

PID

2TI k1 (2Tw + T )

k1 e −Td s T s + 2ξ 0T0 s + 1 PID 0 < ξ0 ≤ 1

2TI k1 (2Tw + T )

k1 −Td s e s

k1 e −Td s (T1s + 1)(T2 s + 1)

*

*

T1 ≥ T2 2 2 0

*

TI*

TD*

-

-

*

T1 −

T 2

T1 −

-

aTI k1

*

aTI k1

*

T 2

T1 + T2 − T

T1T2 T − T1 + T2 4

2ξ 0T0T2 − T

T0 T − 2ξ 0 4

tab. 4.3 Závislost koeficientů přeregulování κ a parametrů α, β

κ

0

0,05

0,10

0,15

0,20

0,025

0,30

0,35

0,40

0,45

0,50

α

1,282 0,948 0,884 0,832 0,763 0,697 0,669 0,640 0,618 0,599 0,577

β

2,718 1,994 1,720 1,561 1,437 1,337 1,248 1,172 1,104 1,045 0,992

Syntéza regulačních obvodů

53

Postup při výpočtu je velmi jednoduchý. Vybere se požadovaná hodnota přeregulování κ, podle tab. 4.3 se určí koeficienty α, β a z nich se určí hodnota koeficientu a: a=

1 , αT + β Td

kde T je vzorkovací perioda a Td je dopravní zpoždění. Pokud známe hodnotu koeficientu a můžeme pomocí tab. 4.2 dopočítat hodnoty jednotlivých koeficientů regulátoru.

Výpočet koeficientů spojitého regulátoru: Metodou požadovaného modelu navrhněte pro soustavu vhodný spojitý regulátor tak, aby hodnota relativního překmitu dosáhla maximálně 5%. Soustava je dána přenosem: GS ( s ) =

1 ⋅ e −14, 4 (80s + 1)(40 s + 1)

V tab. 4.2 najdeme jaký typ regulátoru musíme použít a rovnice pro výpočet jeho koeficientů. Pomocí tab. 4.3 určíme hodnotu koeficientu β pro daný překmit (hodnotu α pro spojitý regulátor nebudeme potřebovat). β = 1,944 Při řešení v programu matlab začneme tím, že zadáme parametry soustavy a pomocí nich zadáme přenos soustavy Gs. k1=1 %zesílení T1=80 %koeficient T1 T2=40 %koeficient T2 Td=14.4 %dopravní zpoždění Gs=tf(k1,[T1*T2 T1+T2 1],'InputDelay',Td) %zadání p řenosu soustavy Transfer function: 1 exp(-14*s) * -------------------3200 s^2 + 120 s + 1

Dále si vypočítáme pomocnou proměnou a. b=1.944 A=1/(b*Td) A = 0.0357

Nyní už známe všechny hodnoty pro výpočet jednotlivých parametrů regulátoru podle rovnic z tab. 4.2. I=T1+T2 I = 120 D=(T1*T2)/(T1+T2) D = 26.6667 P=(A*I)/k1 P = 4.2867

Syntéza regulačních obvodů

54

Pomocí vypočítaných parametrů zadáme přenos regulátoru: Gr=tf([P*I*D P*I P],[I 0]) Transfer function: 1.372e004 s^2 + 514.4 s + 4.287 ------------------------------120 s

Teď již zbývá jen propojit regulátor a soustavu. Jelikož však máme přenos soustavy zadán ve tvaru s dopravním zpožděním budeme jej nejdříve muset aproximovat na tvar bez dopravního zpoždění. K tomu nám souží příkaz pade. Gs1=pade(Gs,3) Transfer function: -s^3 + 0.8333 s^2 - 0.2894 s + 0.04019 -------------------------------------------------------------3200 s^5 + 2787 s^4 + 1027 s^3 + 164.2 s^2 + 5.112 s + 0.04019

Pro porovnání si vykreslíme přechodové charakteristiky původního i aproximovaného přenosu soustavy: step(Gs,Gs1) Step Response From: U(1) 1.2

1

To: Y(1)

Amplitude

0.8

0.6

0.4

0.2

0

-0.2 0

50

100

150

200

250

300

350

400

450

500

Time (sec.)

obr. 4.3 Porovnání přechodových charakteristik soustavy před a po aproximaci

Teď již můžeme regulátor i soustavu propojit do uzavřeného regulovaného obvodu a vykreslit jeho přechodovou charakteristiku (obr. 4.4). Go=Gs1*Gr Gro=feedback(Go,1)

Syntéza regulačních obvodů

55 Step Response From: U(1)

1.2

1

To: Y(1)

Amplitude

0.8

0.6

0.4

0.2

0

-0.2 0

50

100

150

200

250

300

350

400

450

500

Time (sec.)

obr. 4.4 Přechodová charakteristika regulované soustavy

Výpočet koeficientů diskrétního regulátoru: Metodou požadovaného modelu navrhněte pro soustavu vhodný diskrétní regulátor tak, aby hodnota překmitu dosáhla maximálně 5% a vzorkovací perioda T = 5. Soustava je dána stejně jako u výpočtu spojitého regulátoru přenosem: GS ( s ) =

1 ⋅ e −14, 4 (80s + 1)(40 s + 1)

V tab. 4.2 najdeme jaký typ regulátoru musíme použít a rovnice pro výpočet jeho koeficientů v tomto případě půjde o PID regulátor. Při zadání diskrétní formy PID regulátoru si musíme uvědomit, že se zadává v odlišném tvaru než spojitý regulátor:  T z −1 T z . G R ( z ) = k p ⋅ 1 + ⋅ + D⋅ z   TI z − 1 T Pomocí tab. 4.3 určíme hodnotu koeficientů α a β pro daný překmit. α = 0,984; β = 1,944 Při řešení v programu Matlab budeme postupovat stejně jako při výpočtu spojitého regulátoru začneme tím, že zadáme parametry soustavy a pomocí nich zadáme přenos soustavy GS. k1=1 T1=80 T2=40 Td=14.4 T=5

%zesílení %koeficient T1 %koeficient T2 %dopravní zpoždění

Gs=tf(k1,[T1*T2 T1+T2 1],'InputDelay',Td) %zadání p řenosu soustavy Transfer function: 1 exp(-14*s) * -------------------3200 s^2 + 120 s + 1

Syntéza regulačních obvodů

56

Dále zadáme koeficienty α a β a vypočítáme pomocnou proměnou a. a=0.984 %koeficient alfa a = 0.9840 b=1.944 %koeficient beta b = 1.9440 A=1/(b*Td) %promena a A = 0.0357

Nyní už známe všechny hodnoty pro výpočet jednotlivých parametrů regulátoru podle rovnic z tab. 4.2. » TI=T1+T2-T TI = 115 » TD=(T1*T2)/(T1+T2)-(T/4) TD = 25.4167 » kp=(A*TI)/k1 kp = 3.4940

Pomocí vypočítaných parametrů zadáme přenos regulátoru tak, že nejdříve zadáme přenosy jednotlivých složek a ty spolu sečteme. >> I=(T/TI)*tf([-1 0],[-1 1], T) Transfer function: 0.04348 z --------z - 1 Sampling time: 5 >> D=(TD/T)*tf([-1 1],[-1 0], T) Transfer function: 5.083 z - 5.083 --------------z Sampling time: 5 >> Gr=kp*(1+I+D) Transfer function: 21.41 z^2 - 39.02 z + 17.76 --------------------------z^2 - z Sampling time: 5

Jelikož máme přenos soustavy zadán ve tvaru s dopravním zpožděním a pro další práci budeme potřebovat přenos bez dopravního zpoždění je třeba jej aproximovat. K tomu nám souží příkaz pade. » Gs1=pade(Gs,3) Transfer function: -s^3 + 0.8333 s^2 - 0.2894 s + 0.04019 -------------------------------------------------------------3200 s^5 + 2787 s^4 + 1027 s^3 + 164.2 s^2 + 5.112 s + 0.04019

Syntéza regulačních obvodů

57

Teď již zbývá jen propojit regulátor a soustavu. Předtím však ještě musíme spojitý přenos soustavy převést na diskrétní tvar s odpovídající vzorkovací periodou. Pro převod použijeme příkaz c2d. » Gsd=c2d(Gs1,T) Transfer function: 0.0001612 z^4 - 0.0005604 z^3 + 0.0007305 z^2 + 0.00425 z + 0.0004643 --------------------------------------------------------------------z^5 - 2.214 z^4 + 1.659 z^3 - 0.552 z^2 + 0.1245 z - 0.01285

Sampling time: 5

Step Response From: U(1) 1.2

1

To: Y(1)

Amplitude

0.8

0.6

0.4

0.2

0

-0.2 0

50

100

150

200

250

300

350

400

450

Time (sec.)

obr. 4.5 Porovnání přechodové charakteristiky soustavy v diskrétním a spojitém tvaru

Teď již můžeme regulátor i soustavu propojit do uzavřeného regulovaného obvodu a vykreslit jeho přechodovou charakteristiku (obr. 4.6). Go=Gsd*Gr Gro=feedback(Go,1) step(Gro) Step Response From: U(1) 1.2

1

To: Y(1)

Amplitude

0.8

0.6

0.4

0.2

0

-0.2 0

50

100

150

200

250

300

350

400

Time (sec.)

obr. 4.6 Přechodová charakteristika regulačního obvodu

450

Syntéza regulačních obvodů

58

4.3 Metody výpočtu koeficientů diskrétních regulátorů Mezi metody přímého výpočtu diskrétních regulátorů patří: o metody výpočtu koeficientů diskrétních regulátorů na základě kvadratických kritérií (diskrétní formy PID, a všeobecné diskrétní regulátory) o metody ukončení regulačního pochodu za konečný počet kroků – deadbeat (DB– regulátory) o metody založené na algebraické polynomiální teorii – moderní formy výpočtu koeficientů všeobecných diskrétních regulátorů na základě řešení diofantických rovnic. Metody určení koeficientů diskrétního regulátoru vycházející z diskrétní přenosové funkce umožňují dosáhnout lepší kvalitu řízení než metody určení koeficientů regulátoru přepočtem ze spojité na diskrétní formu. Na rozdíl od diskrétních forem spojitých regulátorů může mít všeobecný diskrétní regulátor úplně rozdílnou strukturu už také vzhledem k tomu, že polynomy Q(z) a P(z) jsou vyššího řádu než 3. Přenosová funkce všeobecného diskrétního regulátoru v z-oblasti je ve tvaru: q 0 + q1 z −1 + q 2 z −2 + ... + q n z − n Q( z ) U ( z ) = = GR ( z) = . −1 − 21 −n P( z ) E ( z ) p 0 + p1 z + p 2 z + ... + p n z

(4.1)

Z diskrétní přenosové funkce regulátoru (4.1) vyplývá, že řídící zásah pro všeobecný diskrétní regulátor (pokud p0=1) má tvar u (k ) = − p1u (k − 1) − p 2 u (k − 2) − ... − p m u (k − m) + q 0 e(k ) + q1e(k − 1) + ... + q n e(k − n) . Syntézou diskrétního regulátoru určíme neznámé koeficienty [q0, q1,…, qn; p1, p2,…,pm]. [Kozák, Š. 1999]

Metody ukončení regulačního pochodu za konečný počet kroků Nejběžnějším a kvalitativně novým prvkem v diskrétní teorii automatického řízení je požadavek na ukončení regulačního pochodu za konečný počet kroků: tento způsob řešení se v literatuře označuje jako řízení za konečný počet koků (deadbeat control) a nebo jako konečné časově-optimální řízení a diskrétní přenos uzavřeného regulačního obvodu je racionální lomená funkce. Současné prakticky využívané algoritmy diskrétního řízení používají pro výpočet akčního zásahu ukončení regulačního pochodu klasický přístup, ale také novější přístup založený na tzv. algebraické teorii řízení vzhledem na její větší perspektivy a mohutný numerický aparát řešení tzv. diofantických rovnic, který se při řešení využívá. [Kozák, Š. 1999]

Návrh deadbeat regulátoru v prostředí Matlab Metoda vychází z přenosové funkce řízeného procesu která se získá ze spojitého popisu podle vztahu: GS ( z ) =

z − 1  −1  G ( s )  Z L   z   s 

Syntéza regulačních obvodů

59

(přepočet s tvarovacím členem nultého řádu) V programu Matlab to realizujeme příkazem C2D(num, den, Tv,’metoda přepočtu’). Pokud metodu přepočtu neuvádíme, automaticky je nastavena metoda ZOH tvarovací člen nultého řádu. Přenosová funkce řízeného procesu je potom ve tvaru: b1 z −1 + ... + bM z − M B( z ) = −1 −M A( z ) 1 + a1 z + ... + a M z

G p ( z) =

Pro deadbeat regulátor je přenosová funkce uzavřeného obvodu ve tvaru GY / W ( z ) = p1 z −1 + p 2 z −2 + ... + p M z − M = P( z ) q1=a1q0

p1=b1q0

q2=a2q0

p2=b2q0

………………………… qM=aMq0

protože

∑p

i

=1

∑b q i

0

= 1 odtud

pM=bMq0

q 0 = 1 / ∑ bi = 1 /(b1 +b 2 +... + bM ) = u (0) Charakteristická rovnice diskrétního uzavřeného obvodu je: 1 + G p ( z )⋅G R ( z ) =≡= z M má tedy M pólů v počátku z-roviny, regulovaná veličina za M kroků (M je řád soustavy) dosáhne žádanou hodnotu. Pokud proces řízení obsahuje dopravní zpoždění, potom jsou koeficienty DB regulátoru určeny: q1=a1q0

p1=b1q0

q2=a2q0

p2=b2q0

…………………………… qM=aMq0

pM=bMq0

qM+1=aM+1q0

p1+d=b1+dq0=b1q0

qM+d=aM+dq0=0

pM+d=bM+dq0=bMq0

kde bi jsou koeficienty čitatele pro i=1,2,…,M+d. Protože proces obsahuje dopravní zpoždění počet koeficientů je M+d. Prvních d koeficientů bi bude nulových a od i=1+d jsou koeficienty nenulové. Z nulové hodnoty koeficientů čitatele diskrétního procesu je zřejmé, že nulové budou také koeficienty jmenovatele regulátoru procesu pi, pro i=1,2,…,d a hodnoty koeficientů polynomu jmenovatele A od I=M+1 až po i=M+d. Tedy koeficienty: b1 = b2 = ... = 0 aM +1 = aM + 2 = ... = a M + d = 0 b1+ d = b1

b2+ d = b2

... bM + d =b M

Struktura DB regulátoru představuje speciální typ diskrétního regulátoru ve tvaru: GR ( z ) =

q + q z −1 + q2 z −2 + ... + qn z − n U ( z ) Q( z ) = 0 1 −1 = 1 − P( z ) 1 − p1 z − p2 z −2 − ... − qn z − n E ( z)

Syntéza regulačních obvodů

60

Výstupní regulovaná veličina dosáhne žádanou hodnotu referenční proměnné w(k) za M+d kroků (kde d je dopravní zpoždění, které se z spojitého dopravního zpoždění určí vztahem d=D/T a T je perioda vzorkování) tj. y (k ) = w(k ) = 1 u (k ) = u ( M )

pro k ≥ M + d pro k ≥ M

(ustálená hodnota postoupnosti řídícího zásahu) Charakteristická rovnice uzavřeného regulačního obvodu je: GY / W ( z ) = p1 z −1 + p2 z −2 + ... + p d + M z − ( d + M ) = P( z ) ( p1 = ... = pd = 0) Příklad:

Úkolem je určit koeficienty DB regulátoru přímo v programu Matlab na základě jednoduchých příkazů, zadaný regulovaný proces je spojitý a je popsaný přenosovou funkcí ve tvaru: GS (s) =

B(s) 2s + 1 = e −4 s A( s ) (10 s + 1)(7 s + 1)(3s + 1)

Řešení: Prvním krok je zadání spojitého přenosu procesu a jeho přepočet do Z-oblasti » cit=tf([2 1], 1) Transfer function: 2 s + 1 » men1=tf(1,[10 1]) Transfer function: 1 -------10 s + 1 » men2=tf(1,[7 1]) Transfer function: 1 ------7 s + 1 » men3=tf(1,[3 1]) Transfer function: 1 ------3 s + 1 » Gpd=cit*men1*men2*men3 Transfer function: 2 s + 1 ---------------------------210 s^3 + 121 s^2 + 20 s + 1

Syntéza regulačních obvodů

61

Gpz=C2D(Gpd,4) Transfer function: 0.06525 z^2 + 0.04793 z - 0.007505 -----------------------------------z^3 - 1.499 z^2 + 0.7041 z - 0.09978 Sampling time: 4

Protože proces obsahuje dopravní zpoždění D=4 a perioda vzorkování T=4 v z-oblasti bude dopravní zpoždění d = D/T = 4/4 = 1, teda celá přenosová funkce procesu bude násobená z-1 čím se posunou koeficienty čitatele o 1 řád. Ve výpočtech pro koeficienty regulátoru budou proto nulové koeficienty p(i) a to o tolik jaká je hodnota d. Pro náš případ d = 1, tedy bude nulový první koeficient tj. p = (1 )= 0, ostatní koeficienty budou počítané z posunutých hodnot polynomu B. G pzd

0.06525 z 2 + 0.04793z − 0.007505 −1 = z z 3 − 1.499 z 2 + 0.09978

Výběr čitatele a jmenovatele z přenosu regulovaného procesu » [Bz, Az]=tfdata(Gpz,'v') %vyber citatele a jmenovatele Bz = 0

0.0653

0.0479

-0.0075

1.0000

-1.4986

0.7041

-0.0998

Az =

Upozornění: Matlab indexuje první prvek v matice a nebo ve vektoru od jedné tj. Bz(1) a nezná Bz(0). Koeficienty čitatele regulátoru q0=1/sum(Bz) q0 = 9.4628

Koeficienty jmenovatele regulátoru p1=0 p1 = 0

q1=Az(2)*q0 q1 =

p2=Bz(2)*q0 p2 =

-14.1813 0.6175 q2=Az(3)*q0 q2 =

p3=Bz(3)*q0 p3 =

6.6627 0.4536 q3=Az(4)*q0 q3 = -0.9442

p4=Bz(4)*q0 p4 = -0.0710

Syntéza regulačních obvodů

62

Zadání polynomů regulátoru: Q=[q0, q1, q2, q3] Q = 9.4628 -14.1813

6.6627

-0.9442

P=[p1, p2, p3, p4] P = 0 0.6175

0.4536

-0.0710

Před zadáním přenosu regulátoru můžme provést kontrolu správnosti výpočtů: sum(P) % Kontrola spravnosti vypoctu ans = 1 sum(Q) % Kontrola ustaleneho stavu ans = 1.0000

Přenos DB regulátoru bude mít následující tvar: GR ( z) =

Q( z ) 9.4628 − 14.1813 z −1 + 6.6442 z −2 − 0.9442 z −3 = 1 − P( z ) 1 − 0.6175 z − 2 − 0.4536 z −3 + 0.070 z − 4

Závěr

63

5 Závěr Cílem mojí práce bylo zpracování elektronické příručky pro podporu výuky automatického řízení. Tato příručka obsahuje řešení typových příkladů z oblasti automatizace pomocí programového prostředí Matlab. V první části jsem popsal program Matlab a jeho toolboxy Simulink, Control System, System Identification, Signal Processing a Optimization. Dále jsem zjistil možnosti jejich využití při řešení úkolů z oblasti automatického řízení. Při své další práci jsem pak využíval hlavně základní funkce programu Matlab, toolbox Control System a Simulink. Druhá kapitola je zaměřena na oblast analýzy systémů. Zde jsem vybral příklady na zadávání matematických popisů a systému a převody mezi nimi, dále pak příklady na výpočet přechodových, impulsních a kmitočtových charakteristik a příklady na kontrolu stability systému. Pro vybrané příklady jsem vypracoval řešení za použití programu Matlab. U většiny příkladů je uvedena řešení jak pro spojité tak pro diskrétní systémy. V třetí kapitole zaměřené na oblast syntézy jsou uvedeny postupy pro návrh spojitých a diskrétních regulátorů v prostředí Matlab. Pro výpočet spojitých regulátorů jsem vybral příklady řešené metodou Ziglera-Nicholse a metodou požadovaného modelu. Pro diskrétní regulátory uvádím dva postupy řešení. První popsaný způsob výpočtu diskrétního regulátoru využívá metodu požadovaného modelu pro návrh diskrétní formy PID regulátoru. Druhý postup využívá metodu ukončení regulačního pochodu za konečný počet kroků tzv. Deadbeat regulátor. Posledním úkolem bylo vytvoření elektronické příručky v hypertextové podobě, která bude přístupná na Internetu a měla by sloužit studentům. Jsem přesvědčen, že v této práci jsou splněny splnil všechny body zadání bakalářské práce. Vytvořená elektronická příručka je zcela funkční a může plnit svoji funkci doplňkového studijního materiálu.

Literatura

64

6 Literatura DUŠEK F. 2000 Matlab a Simulink úvod do používání, Univerzita Pardubice, Bratislava, 2000, ISBN 80-7194-273-1. FARANA, R. AJ, 1996. Programová podpora simulace dynamických systémů. Sbírka řešených příkladů. 1 vyd. Ostrava: KAKI 1996, 114 s. ISBN 80-02-01129-5. KOZÁK, Š. 1999. Matlab – Simulink II. 1.vydání, STU Bratislava, 1999. 141 s. ISBN 80-227-1235-3 KUPCZAK, M. 2002. Programová podpora analýzy regulačních obvodů pomocí programu MATLAB. Ostrava: katedra ATŘ-352 VŠB-TU Ostrava, 2002. 71 s. Diplomová práce, vedoucí: Vítečková, M. MIZERA, R. 2002. Programová podpora syntézy regulačních obvodů pomocí programu MATLAB-Simulink. Ostrava: katedra ATŘ-352 VŠB-TU Ostrava, 2002. 59 s. Diplomová práce, vedoucí: Vítečková, M. MODRLAK, O. VOLEJNÍK, O. Stručný manuál Matlabu pro předměty teorie řízení. [online]. Dostupné na URL: (14.4.2003) NOSKIEVIČ, P. 1999. Modelování a identifikace systémů. Ostrava: Montanex 1999. 276 s. SIGMON, K. 1992. MATLAB Primer. University of Florida, 1992. on line na URL: (14.4.2003) SLOVÁK, T. 2002. Využití simulačního programu SIPRO pro identifikaci a prezentaci výsledků v prostředí Internetu. Ostrava: katedra ATŘ-352 VŠB-TUO, 2002. 96 s. Bakalářská práce, vedoucí: Wagnerová, R. SMUTNÝ, L. Interní předpisy pro vypracování informační stránky diplomové/bakalářské práce. Ostrava : VŠB-TU Ostrava, 2001 [cit. 2001-01-22]. Dostupný z www na URL: ŠVARC, I. 1992. Teorie automatického řízení I. Brno: FS VUT, 1992. 210 s. ŠVARC, I. 1993. Teorie automatického řízení II. Brno: FS VUT, 1992. 231 s. THE MATHWORKS, INC. 1996. Control System Toolbox User’s Guide. Second printing for MATLAB 5. THE MATHWORKS, INC. 1996. Optimization Toolbox User’s Guide. Second printing for MATLAB 5. THE MATHWORKS, INC. 1996. Signal Processing Toolbox User’s Guide. First printing for MATLAB 5. THE MATHWORKS, INC.2001. Systém Identification Toolbox User’s Guide. Fifth printing. THE MATHWORKS, INC.1999. Using Simulink. Revised for Simulink 3 (Relase 11). VÍTEČKOVÁ, M. Seřízení regulátorů metodou inverze dynamiky 1. vyd. Ostrava : VŠB-TU Ostrava, 1998. 56 s. ISBN 80-7078-628-0 WAGNEROVÁ, R. – MINÁŘ, M. AJ. 2000. Syntéza regulačních obvodů. Ostrava: VŠBTUO, KATŘ, 2000. Dostupné na URL: (14.4.2003) WAGNEROVÁ, R. – MINÁR, K. AJ. 2000. Prezentační a výukový modul pro oblast analýzy regulačních obvodů v prostředí INTRANETU. Ostrava: VŠB-TUO, KATŘ, 2000.dostupné na URL: (14.4.2003)