Szabályozástechnika - Matlab alapismeretek
A VIK Wikiből
A mechanikai lengőrendszer leírása
clear all; % Az összes változó törlése a munkatérből
clc % A Command Window kiürítése
% A rendszer paraméterei:
m=0.5; % Test tömege [kg]
k=0.2; % Rugóállandó [N/m]
b=0.5; % Csillapítás [Ns/m]
% A rendszer differenciálegyenlete:
% F = mx'' + bx' +kx
% Ugyanez s tartományban:
% F(s) = m*s^2*x(s) + b*s*x(s) + k*x(s)
% A rendszer átviteli függvénye:
% W(s) = x(s)/F(s) = 1/(ms^2+bs+k)
% Állapotteres leírás mátrixokkal
% Állapotváltozók, be és kimenet:
% x_ = [x x']'
% u = F
% y = x
%
% Az állapotegyenletek:
% x' = x'
% x'' = -(k/m)x - (b/m)x' + (1/m)F
%
% Az állapotteres leírás:
% x_' = Ax_ + Bu
% y = Cx_ + Du
Állapotteres reprezentáció
% Az állapotteres leírás mátrixai:
A = [0 1; -k/m -b/m]
B = [0; 1/m] % Másképpen: B=[0 1/m]'
C = [1 0]
D = 0
% x hatása x''-re
A(2,1) % Az A mátrix 2. sorának 1. eleme (index: sor, oszlop)
% az x''-re vonatkozo egyenlet
A(2,:) % Az A mátrix 2. sora
% Speciális mátrixok
E=eye(4) % 4x4-es egységmátrix
O=ones(2,4) % 2x4-es, 1 elemeket tartalmazó mátrix
Z=zeros(4,2) % 4x2-es nullmátrix
% Mátrixok mérete
size(A) % A mátrix mérete: [sorok száma, oszlopok száma]
size(B,1) % B vektor sorainak száma
size(C,2) % C vektor oszlopainak száma
length(B) % B vektor hossza
% Mátrixfüggvények
eig(A) % A mátrix sajátértékei
inv(A) % A mátrix inverze
sys=ss(A,B,C,D) % SS objektum létrehozása a mátrixokból
fieldnames(sys) % Az SS objektum tulajdonságai
A=sys.a % A sys objektum a tulajdonsága (A mátrix)
B=sys.b % A sys objektum b tulajdonsága (B mátrix)
C=sys.c % A sys objektum c tulajdonsága (C mátrix)
D=sys.d % A sys objektum d tulajdonsága (D mátrix)
[A,B,C,D]=ssdata(sys) % Az állapotteres leírás mátrixainak kinyerése az objektumból
Átviteli függvény reprezentáció
% Átviteli függvények megadása polinomokkal
num=1 % A számláló egy skalár
% A nevező egy polinom: ms^2+bs+k
% Polinomokat vektorként írunk le, aminek elemei az együtthatók a hatványok csökkenő sorrendjében
den=[m b k]
polyval(den,1) % s=1 behelyettesítése a den polinomba
polyder(den) % A den polinom s szerinti deriváltja
p=roots(den) % A den polinom gyökei (vektorként)
% Polinom összeállítása gyökökből
p1=[1 -p(1)] % p1=s-p(1)
p2=[1 -p(2)] % p2=s-p(2)
conv(p1,p2) % Polinomszorzás: conv() függvény
% Három polinom szorzata: conv(conv(p1,p2),p3)
poly(p) % Polinom létrehozása a gyökök vektorából
w=tf(num,den) % TF objektum létrehozása a számláló és nevező polinomjából
fieldnames(w) % A TF objektum tulajdonságai
num2=w.num{1} % Számláló és nevező polinom: num és den mezők
den2=w.den{1} % Figyelem! A {1} kötelező!
[num2,den2]=tfdata(w) % Számláló és nevező kinyerése
% Megadás szimbolikus 's' segítségével
s=tf('s') % A bal oldali 's' egy szimbolikus változó
wsym=1/(m*s^2+b*s+k)
Zérus-pólus-erősítés reprezentáció
% Pólusok és zérusok
p=roots(den) % p vektor az előbbi nevező gyökeit tartalmazza
p1=p(1) % p1 az első pólus
abs(p1) % p1 abszolút értéke
phase(p1) % p1 fázisszöge
real(p1) % p1 valós része
imag(p1) % p1 képzetes része
z=roots(num) % z vektor a számláló gyökeit tartalmazza
% Mivel a zpk alakban a legmagasabb fokszámú együttható mind a számlálóban,
% mind a nevezőben 1, így a k erősítés értéke
k=num(1)/den(1)
% Zérus-pólus-erősítés reprezentáció
syszp=zpk(z,p,k) % ZPK objektum létrehozása
fieldnames(syszp) % A ZPK objektum mezői
syszp.z{1} % Zérusok és pólusok
syszp.p{1} % A {1} kötelező!
syszp.k % Erősítés - nem kell {1}!
Áttérés a reprezentációk között
% Bármely LTI objektum (SS,TF,ZPK) konstruktora meghívható úgy is, hogy paramétere egy másik LTI objektum
tf(sys) % Állapotteres leírás -> átviteli függvény
zpk(sys) % Állapotteres leírás -> zérus-pólus-erősítés
ss(w) % Átviteli függvény -> állapotteres leírás
zpk(w) % Átviteli függvény -> zérus-pólus-erősítés
tf(syszp) % Zérus-pólus-erősítés -> átviteli függvény
ss(syszp) % Zérus-pólus-erősítés -> állapotteres leírás
Vizsgálati módszerek az időtartományban
% Klasszikus vizsgálójelek
step(sys) % Ugrásválasz
step(sys,w,syszp) % Több rendszer válaszának ábrázolása
legend('SS','TF','ZKP') % Jelmagyarázat
impulse(w) % Impulzusválasz
% Jobb klikk --> Characteristics menüpont:
% Stady State - A függvény végértéke
% Rise Time - A függvény felfutási ideje
% Settling Time - A függvény 2%-os beállási ideje
% Peak Response - A függvény maximuma és túllövése
x0=[2,0] % Kezdeti állapot: x0=2m, v0=0m/s
initial(sys,x0) % A rendszer válasza 0 bemenet mellett
% Válasz tetszőleges jelre
t=[0:0.01:10]; % Idővektor: 0,0.01,0.02,...9.99,10
u=4*sin(5*t); % u vektor: 4(sin 0), 4(sin 0.05), ...
lsim(w,u,t) % Az u(t) gerjesztőjel és az arra adott válasz a t=0, t=0.01, t=0.02... időpontokban
% Vizsgálati módszerek a frekvenciatartományban
pzmap(syszp) % Pólus-zérus elrendezés
bode(w) % Bode-diagram
nyquist(w) % Nyquist-diagram
Szabályozási kör összeállítása
% A lengőrendszer átviteli függvénye
wdamp=tf(1,[m b k])
% Legyen az erőt szolgáltató beavatkozó szerv egy egytárolós tag, melynek időállandója
Tact=1
% A beavatkozó szerv átviteli függvénye
wact=tf(1,[Tact 1])
% A szakasz átviteli függvénye a két rendszer soros kapcsolása
wp=wdamp*wact % Másképpen: wp=series(wdamp,wact)
% Vizsgáljuk meg a szakaszt!
% A szakasz pólusai
roots(wp.den{1})
% A damp() utasítás visszaadja nemcsak a szakasz pólusait, de a pólusokhoz tartozó
% csillapítatlan sajátfrekvencia (w0) és csillapítás (xi) értéket is
damp(wp)
% A szakasz ugrásválasza
step(wp)
% Egyszerű P-szabályozó méretezése túllövésre
% Alkalmazzunk P szabályozót (a szabályozó egy egyszerű erősítés), majd
% a zárt körben állítsunk be ~10%-os túllövést!
% Ehhez először rajzoltassuk ki a gyökhelygörbét, és keressük meg azt az
% erősítést, amely mellett a túllövés (overshoot) körülbelül 10% lesz!
rlocus(wp)
% Az erősítés értéke a gyökhelygörbéről leolvasva:
Ap=0.0696
% Így a szabályozó átviteli függvénye wc(s) = Ap:
wcp=tf(Ap,1)
% A felnyitott kör átviteli függvénye w0 = wc*wp
w0p=wcp*wp
% Ellenőrizzük a fázis-, és erősítéstartalékot a margin függvény segítségével
margin(w0p)
% Állítsuk össze a zárt kört és vizsgáljuk meg a tulajdonságait!
% A zárt kört a feedback utasítás használatával állíthatjuk elő. Ennek első paramétere az
% előre vezető ág (LTI rendszer, itt átviteli függvény), második paramétere a visszacsatoló
% ág (LTI rendszer, itt átviteli függvény), harmadik paramétere pedig a visszacsatolás előjelét
% határozza meg. Mivel egységnyi merev visszacsatolást alkalmazunk, ezért a visszacsatoló ág
% átviteli függvénye w(s)=1, azaz tf(1,1), a visszacsatolás előjele pedig negatív, ezért a harmadik
% paraméter egy negatív szám (szokásos a -1 választás).
wclp=feedback(w0p,tf(1,1),-1)
% Egységnyi merev, negatív visszacsatolás esetén egyszerűbben is írható: wclp=feedback(w0p,1)
% Vizsgáljuk meg a zárt kör pólusait mind numerikusan, mind grafikusan!
damp(wclp)
pzmap(wclp)
% Rajzoljuk ki a zárt kör ugrásválaszát!
step(wclp)
% A zárt kör ugrásválaszának értékét numerikusan is meghatározhatjuk a Wcl(s=0) helyettesítéssel
% vagy a dcgain függvénnyel (vigyázat, csak nulladik típusú kör esetén működik!)
K=polyval(wclp.num{1},0)/polyval(wclp.den{1},0)
K=dcgain(wclp)
% Határozzuk meg a beavatkozó jelet zárt körben! Ehhez állítsuk elő a zárt körbeli r->u átviteli
% függvényt, amihez szintén a feedback utasítást használjuk. Ezúttal az előre vezető ágban wc,
% míg a visszacsatoló ágban wp helyezkedik el, az előjel pedig továbbra is negatív.
wu=feedback(wcp,wp,-1)
step(wu)
