Szabályozástechnika - Matlab alapismeretek

← Vissza az előző oldalra – Szabályozástechnika

 

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)