Szabályozástechnika - Laborjegyzetek
HalacS által feltöltött jegyzet3.txt (3. labor jegyzet (utf8)), jegyzet4.txt (4. labor jegyzet (utf8)), jegyzet5.txt (5. labor jegyzet (utf8)) fájlokból inkább Wiki-aloldalt készítettem, így könnyen szerkeszthető-bővíthető.--Haraszin Péter (vita) 2013. május 16., 09:50 (UTC)
Tartalomjegyzék
3. labor
1) stabilitás vizsgálat a ZÁRT rendszer pólusai alapján
Azért zárt, mert a ki és bemeneten kívül mást nem látunk belőle.
(s-z1)(s-z2)...
H(s)=k * --------------------
(s-p1)(s-p2)...
A zárt rendszer pólusainak valós része < 0 <=> a rendszer stabil
T(s)= s+5/(s^5+3s^4+4s^3+10s^2+5s-10)
Matlab:
num = [1 5]; den = [1 -3 4 10 5 -10]; T = tf(num,den);
Egy rendszer pólusai, a nevező gyökei ezért: Matlab:
p = roots(den) step(T)
Ebből érezhető, hogy nem stabil, de ennyi indoklás nem elég, mert lehet hogy később (amit nem ábrázol) már stabillá válik.
2) Nyquist stabilitási kritérium
Itt már van egy negatív visszacsatolás is a rendszerben, nem csak magára van hagyva a francba a rendszer. A bemenetet a kimenet függvényében húzza el a bemenetet, hogy azt kapjuk a kimeneten amit szeretnénk. Azaz a bemenetre a hiba kerül. Ha a hiba 0, nem csinál semmit.
Felnyitott kör: a visszacsatolást elhadjuk. Zárt szabályozási kör: visszacsatolással.
Tuy = L/(-1(L*(-1))) = L / (1+L)
Ha csak egy szimpla negatív visszacsatolás van és egy doboz fölül, aminek ismerjük az átvitelét, akkor jó a nyquist kritérium: A ZÁRT rendszer stabil (T), ha a NYITOTT KÖRNEK (L), PONTOSAN annyi labilis pólsa van ahányszor a felnyitott kör nyquist diagrammja megkerüli a -1+0j pontot az óra járásával ellentétes irányban.
Labilis pólus: a valós része pozitív
felvágja egy diagrammra a pólusokat meg a zérusokat: pzmap(T) pólus: x zérus: o (az ikszecskék a bal félsíkon legyenek akkor stabil)
1. példa
-5
L(s) = --------------- <<-- a felnyitott kör átvitele
(1+10s)(1+s)
L
T = ----- <<-- zárt kör átvitele
1+L
Matlab:
s = zpk('s')
L=(-5)/((1+10*s)*(1+s))
nyquist(L)
- --> pontosan 0 megkerülés 0 labilis pólussal -> a zárt rendszer stabil
Megnézzük azért az egész rendszerre is a stabilitást:
Matlab:
T=L/(1+L) pzmap(T)
- --> a jobb félsíkot már nem is ábrázolta, mert minden ikszecske a bal oldalt van -> oké
A 2. ábrán gyorsabban áll be stabilba. A különböző erősítés nem érdekes, mert arra lehet konstans erősítést rakni azt jóvan.
Matlab
figure(1), step(L), figure(2), step(T)
-5
L(s) = ------------------
(1-10s)(1+0,1s)
1-10s=0 P1 = 1/10 > 0 => instabil 1+0,1s=0 s=-10
Matlab:
L2=-5/((1-10*s)*(1+0.1*s)) nyquist(L2)
- A*sin(wT) => G(w)*A*sin(wt+e(w))
- A szög megmondja mennyire tér el a sin argumentuma a távolság pedig az erősítést.
Matlab:
T2=L2/(1+L2) T2=minreal(T2) figure(1), step(L2), figure(2), step(T2)
- instabil rendszerből egy negatív visszacsatolással stabil rendszert csináltunk.
A stabilitás jó hogy van, de kell valami mérőszám is hozzá.
- fázistartalék
- erősítési tartalék
- modulus tartalék
- késleltetési tartalék
L most stabil. Az a jó Nyquist-diagramm, ami a -1+0j és a 0 között megy el. A rendszer legyen nagyon stabil, azaz ne mennyen a nq diagramm túl közel a -1+0j ponthoz. Ez az erősítési tartalék (gm = gain margin) gm = 1/G(w) (de lehet máshogy is definiálni, csak ő ezt szereti)
A hibajelet többszörözve is belehet vinni a bemenetre (összegző és bemenet közé erősítő), ekkor gyorsul a stabilba állás (nagyobb hibát hazudunk be neki), de ezzel isntabilba is ellehet vinni a rendszert.
1
L = ---------------
(0.5+s)(s^2+2s+1)
Matlab:
L = 1/((0.5+s)*(s^2+2*s+1)) nyquist(L)
- --> wpi = 1,41
- Re = -0,225
- Im = 0
- -> g=1/0,225=4,444
Matlab:
TT1=4*L/(1+4*L) pzmap(TT1) % --> csak a negatív félsíkon vannak pólusai (zérusai??) a rendszernek TT2=5*L/(1+5*L) pzmap(TT1) % --> itt már a jobb félsíkon is vannak pólusok
Fázisban is érdekes, hogy milyen messze vagyunk a -1+0j ponttól. Lehet olyan perverz is a rendszer, hogy erősítési tartalékunk sok van, de fázis tartalék kevés. Ez az a szám, ami a legközelebbi pontnak a vízszinteshez viszonított eltérése fokban.
Matlab:
[gm, pm, wg, wc] = margin(L) % ^ ^ ^ ^ % | | | | % | | | vágási kör frekvencia % | | wpi % | phase margin (a 60 fok már egész jó) % gain margin
Bode diagrammból:
erősítési tartalék: ahol az alsó lemegy 180 fokra, azon a ponton a felső diagrammból 1/G(w) fázis tartalék: ahol a felső metszi az x tengelyt, ott az alsó 180-fi(w)
Modulus tartalék: a -1+0j pontból az első egység kör ami elmetsze a Nyquist diagrammot. (az erősítést és a fázisrontást kombinálva elszaródás kb) számolása:
- ro+ = min|1+L(iw)|
- w
késleltetési tartalék: ennyi holtidőt lehet plusszba hozzáadni a rendszerhez hogy még stabil maradjon
Matlab:
[mag, phase, w] = bode(L) mag = mag(:) phase = phase(:) Tbl = [mag, phase, w]
- ->lett egy táblázat, aminek az oszlopai: G(w), E(w), w
- ebből erősítési tartalék keresése: addig görgetjük a táblázatot, még meg nem találjuk, hogy G(w), -180, wpi (itt konkrétan: 0.2503 -176.0386 1.3394). Ebből tudjuk, hogy az erősítési tartalék G(w)=0.25 -> 1/0.25 = 4
fázistartalék: addig keresünk az előbbi táblázatban, amég olyat nem találunk, hogy 1, pw, wc (itt konkrétan: 0.9927 -108.3126 0.5715) Ebből a fázistartalék: pw=-108 -> -180-(-108)=72
Modulus tartalék keresése:
Matlab:
[mag2, phase2, w2] = bode(1+L) mag2 = mag2(:) ro=min(mag2) % --> ro = 0.6317
4. labor
2010-10-20 - 4. Labor
Most stabil folyamatokhoz (P) akarunk szabályózt (C) tervezni, úgy hogy a P gyorsabb meg ilyesmibb legyen. A PID szabályzó elég a nem kritikus dolgokhoz: szoba hőmérséklet például) Ez a rész a gyakban használható.
L = C*Pc=?
PID = Proportional (arányos tag), D.. (deriváló tag), I.. (integráló tag)
P -> (5) => 5*e (e=hiba)
D -> a hibajel tendenciájának megfelelően "előredolgozik", így ha hirtelen váltásvan a hibajelben, akkor gyorsabban reagál
I -> kiintegrálja a hibajelet. Megnézi a múltban mennyi votl a hiba. Megakadájozza, hogy offsettel álljon be a kimenet. Ez azét van mert egy idő után nem reagálna a rendszer az akkorra már túl kis arányos tag miatt.
P (arányos) tag csinálás:
P=k (k konstans)
D (deriváló) tag csinálás:
Y(s) = H(s)*U(s)
D(s)=(s-Kd)/1 <-- ilyet nem tudunk csinálni, mert a számlálónak legalább annyinak kéne lennie mint a nevezőnek
s*Kd
D(s)=-------- <-- közelítő deriváló tag
1+s*Tn
Matlab:
s=tf('s')
C=1+s*10
Transfer function: 10 s + 1
Matlab:
step(C) % hiba... ??? Error using ==> rfinputs % Not supported for non-proper models.
s*Ks 1+s*Tn + n*K1 1+s*Tn
PD(s) = 1 + ---------- = ------------------ ~ ----------------
1+s*Tn 1+s*Tn 1+s* (Td/n)
n1=3/2 n2=5
Matlab:
C1= (1+3*s)/(1+2*s) Transfer function: 3 s + 1 ------- 2 s + 1 C2=(1+s*10)/(1+2*s) Transfer function: 10 s + 1 -------- 2 s + 1 figure(1), step(C1, 'b', C2, 'r') figure(2), bode(C1, 'b', C2, 'r')
Az n paraméterrel tetszőlegesen választható meg a pólus.
Integráló szabályozó:
1+s*Ti
PI: 1+1/s*k = 1+1/(s*Ti) = --------
s*Ti
k: mennyire büntetjük a hibát
Matlab:
C3=1+1/(s*2) Transfer function: 2 s + 1 ------- 2 s step(C3)
Együttes PID szabályozó:
1+s*Ti 1+s*Td
PID = k * ---------- * ---------------
s*Ti 1+s*(Ts/m)
Matlab:
P=1/((1+10*s)*(1+s)*(1+0.2*s))
Transfer function:
1
-----------------------------
2 s^3 + 12.2 s^2 + 11.2 s + 1
step(P)
kb 50mp alatt áll be oda ahova szeretne ez a P folyamat. p1 = -1 p2=-5
Az, hogy milyen gyorsan áll be, attól függ, hogy milyen messze vannak a pólusok az origótól. Ha közel van, lassabban áll be. e^(p^i*t) -> e^(-t)=1/e^(t) Azt szeretnénk, ha a pólusok minél nagyobbak és negatívak lennének.
Megválasztható paraméterek a PID szabályzóban: k, Ti, Td, n Ti és Td (integráló és deriváló időállandók) megválasztásával kiüthetjük a nem tetsző, lassító pólusokat.
Azt szeretjük, ha a nevező 1+T*s alakú, mert ebből leolvasható az időállandó. Akkor van baj, ha a T számok nagyok, mert akkor lesz a rendszer lassú.
Matlab:
n = 1; k = 1; Ti=10; Td = 1; C=k*(1+s*Ti)/(s*Ti)*(1+s*Td)/(1+s*Td/n) Transfer function: 10 s^2 + 11 s + 1 ----------------- 10 s^2 + 10 s
L=C*P
Transfer function:
10 s^2 + 11 s + 1 ------------------------------------------- 20 s^5 + 142 s^4 + 234 s^3 + 122 s^2 + 10 s
L=zpk(L)
Zero/pole/gain:
0.5 (s+1) (s+0.1)
-----------------------
s (s+5) (s+1)^2 (s+0.1)
L=minreal(L) % leegyszerűsíti a törtet
Zero/pole/gain:
0.5
-------------
s (s+5) (s+1)
T=L/(1+L) % negatívan visszacsatoljuk
Zero/pole/gain:
0.5 s (s+5) (s+1)
---------------------------------------------
s (s+5) (s+5.025) (s+1) (s+0.8595) (s+0.1158)
figure(2), step(T)
n =
2
C=k*(1+s*Ti).... % újra számoljuk
Transfer function: 10 s^2 + 11 s + 1 ----------------- 5 s^2 + 10 s L=C*P % újra számoljuk
Transfer function:
10 s^2 + 11 s + 1
------------------------------------------
10 s^5 + 81 s^4 + 178 s^3 + 117 s^2 + 10 s
T=L/(1+L) % újra számoljuk
Transfer function: 100 s^7 + 920 s^6 + 2681 s^5 + 3209 s^4 + 1565 s^3 + 227 s^2 + 10 s ----------------------------------------------------------------------------------------------------------- 100 s^10 + 1620 s^9 + 10121 s^8 + 31276 s^7 + 51758 s^6 + 45953 s^5 + 20458 s^4 + 3905 s^3 + 327 s^2 + 10 s
[gr,pm,wg,wc]=margin(L)
gr =
70.0000 % erősítési tartalék --> ha 70x-re erősítjük, akkor is stabil marad (tip.nem kell)
pm = % fázistartalék (ez lesz érdekes)
85.9980
wg =
3.1623
wc =
0.0998
[mag, phase, w] = bode(L); g = margin(mag, phase-60, w) % ha csak egy paramétert adunk meg visszatérésnek a marginnak, akkor csak az elsőt adja vissza
g = % ez az erősítési tartalék!
8.3651
k=g;
C=k*(1+s*Ti)... % újra számoljuk
Transfer function:
83.65 s^2 + 92.02 s + 8.365
---------------------------
5 s^2 + 10 s
L=C*P % újra számoljuk
Transfer function:
83.65 s^2 + 92.02 s + 8.365
------------------------------------------
10 s^5 + 81 s^4 + 178 s^3 + 117 s^2 + 10 s
T=L/(1+L) % újra számoljuk a negatív visszacsatolást
Transfer function: 836.5 s^7 + 7696 s^6 + 2.243e004 s^5 + 2.684e004 s^4 + 1.309e004 s^3 + 1899 s^2 + 83.65 s ------------------------------------------------------------------------------------------- 100 s^10 + 1620 s^9 + 10121 s^8 + 3.201e004 s^7 + 5.853e004 s^6 + 6.57e004 s^5 + 4.409e004 s^4 + 1.543e004 s^3 + 1999 s^2 + 83.65 s
step(T)
Ez már sokkal gyorsabban beáll.
A k-t minél nagyobbra állítjuk annál gyorsabb lesz a rendszer, de ezzel együtt a beavatkozó jelek is megnőnek, ami nem jó ha korlátozni kell a beavatkozó jeleket. Azt fogjuk csinálni az egész félévben, hogy a k-t úgy választjuk meg, hogy maradjon meg a 60 foknyi fázistartalék, mert az nekünk valamiért marha jó lesz.
A tervezésnek mindig az a szempontja, hogy a legnagyobb időállandót mindig az integráló tagban kell elsütni. Ebből: Ti = a legnagyobb T, Td = a második T, n = 2..15 (ha n nagy, nagy a beavatkozó jel, ha n kicsi, lassabb rendszer)
Ha nem kell deriváló tag (PI szabályozó kell), akkor eldobjuk a szabályozó képletéből a második törtet, azt jóvan. Ha nem kell integráló tag (PD szabályozó kell), akkor fordítva hadjuk el, viszont ilyenkor a második legnagyobb T lesz a Td (ez amúgy ugye a Ti lenne ha lenne integráló tag)
(A PID szabályózóval lehet pénzt is keresni...)
Állapotegyenletesdi:
x'= A*x + B*u y = C*x + D*u
Matlab:
P=6/((s+1)*(s+2)*(s+3))
Transfer function:
6
----------------------
s^3 + 6 s^2 + 11 s + 6
P=ss(P) % a P átviteli fv-t áttranszformálja állapot mátrixokra
a =
x1 x2 x3
x1 -6 -2.75 -1.5
x2 4 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 0 1.5
d =
u1
y1 0
Continuous-time model.
Olyan átviteli függvény kéne, ami a köv. pólusokat tudja:
Pk = [-6 ; -4+i*4 ; -4-i*4]
Pk = -6.0000 -4.0000 + 4.0000i -4.0000 - 4.0000i
k=acker(P.a, P.b, Pk)
k =
8.0000 17.2500 46.5000 % ezzel a vektorral elvileg az előbbi pólusokat kapnánk
x' = A*x + B(u-K^T*x) alfa=(a-B*K^T)*x + B*u
Matlab:
A1 = P.a-P.b*k
A1 =
-14 -20 -48
4 0 0
0 1 0
B1 = P.b
B1 =
1
0
0
C1 = P.c
C1 =
0 0 1.5000
D1 = P.d
D1 =
0
P2 = ss(A1,B1,C1,D1)
a =
x1 x2 x3
x1 -14 -20 -48
x2 4 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 0 1.5
d =
u1
y1 0
Continuous-time model.
pzmap(P2)
Az új rendszer az előzőhöz képest torzít:
lp=dcgain(P) % megnézi mennyi a statikus átvitele a fv-nek (gyak az egységugrásra a váalsz)
lp =
1
lp2=dcgain(P2)
lp2 =
0.0313
Kell egy kis korrekció, hogy a statikus erősítés ugyanaz legyen, mint volt:
P3=lp/lp2*P2
a =
x1 x2 x3
x1 -14 -20 -48
x2 4 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 0 48
d =
u1
y1 0
Continuous-time model.
step(P,'r',P2,'b',P3,'g')
5. labor
2010.11.03
Nulladrendű tartószerv: két mintavételezés között konstansnak vesszük a jelet (a félév során mindenhol nulladrendű tartószervekkel fogunk dolgozni)
z-(Laplace-) transzformáltja egy jelnek: int 0..inf y(t)*e^(-st)dt ahol z:= e^(-s) így már ez sum k=0..inf y[k]*z^(-k) lett.
impulzus: sum k=0..inf y[k]*z^(-k) = y[0]*1 = 1
Y(s) = H(s) * U(s) g(z): mintavételes ekvivalense a rendszernek Y(z)=G(z)*U(z) Y(z)=G(t)*1 ha impulzus gerjesztés van
A mintavételi idő változásával együtt változik majd a rendszer modellje is.
y(0) = lim k->0 y(k*Ts) = lim z->inf k=inf vv(z) y(inf)=lim k=inf y(k*Ts)=lim z->1 (1-z^(-1))*vv(z)
Mathlabban: h(1) ==> c2d(H,Ts, 'zoh') %zoh = zero order...
() () ()
H(s) = --------------------
(s-p1)(s-p2)...
p1~>z1=e^(p1*Ts)
z=e^(p*Ts) = e^((a9i*n)Ts) = e^(a*Ts) * e^(i*b*Ts)
<=1 |..|=1
Stabilitási kritériumok: FI: a pólusok legyenek negatívak DI-stabilitás: a pólusok absz. értéke legyen negatív
e^(-s)
P(s) = ------------------
(1+10s)(1+5s)
z=e^(s*Ts) => y^(-k)
Matlabbal hogyan kell megetetni:
s=tf('s')
% Transfer function:
% s
Pp=1/((1+10*s)*(1+5*s))
% Transfer function:
% 1
% -----------------
% 50 s^2 + 15 s + 1
Gp=c2d(Pp,1) % 1 a mintavételi idő
% Transfer function:
% 0.009056 z + 0.008194
% ----------------------
% z^2 - 1.724 z + 0.7408
%
% Sampling time: 1
Gp=zpk(Gp) % hogy lássuk a zérusokat és pólusokat
% Zero/pole/gain:
% 0.0090559 (z+0.9048)
% ---------------------
% (z-0.9048) (z-0.8187)
%
% Sampling time: 1
z=tf('z',1) % ez egy 1 mintavételi idejű diszkrét átviteli fv
% Transfer function:
% z
%
% Sampling time: 1
G=Gp/z
% Zero/pole/gain:
% 0.0090559 (z+0.9048)
% -----------------------
% z (z-0.9048) (z-0.8187)
%
% Sampling time: 1
step(Pp,Gp) % legyen egy rajzon a két ugrásválasz
step(Pp,Gp,G) % nézzük meg a holtidőt is
- a holtidőt bele tudtuk tenni a folyamatba, úgy hogy majd szépen tudjuk kezelni is
z-e^(-Ts/Ti)
PI: A * ------------------
z-1
z - e^(-Ts/Td) z - e^(-Ts/Td)
PD: A * --------------------- = -----------------
z-e^(u*(-Ts/Td)) z
z-e^(-Ts/Ti) z-e^(-eTs/Td)
PID: A * ------------------ * ------------------
z-1 z
C1=(z-0.9048)*(z-0.8187)/(z-1)/z % a két konstans az előző zpk-ból van
% Transfer function:
% z^2 - 1.724 z + 0.7408
% ----------------------
% z^2 - z
%
% Sampling time: 1
L1=C1*G % a felnyitott kör átvitele
% Zero/pole/gain:
% 0.0090559 (z-0.9048) (z-0.8187) (z+0.9048)
% ------------------------------------------
% z^2 (z-1) (z-0.9048) (z-0.8187)
%
% Sampling time: 1
[mag,phase,w] = bode(L1); % megnézzük mennyi a fázis
k=margin(mag,phase-60,w) % lecsalunk a fázisból 60-at, mert kell nekünk
k =
15.4261
C=k*C1
% Transfer function:
% 15.43 z^2 - 26.59 z + 11.43
% ---------------------------
% z^2 - z
%
% Sampling time: 1
T=(C*G)/(1+C*G) % a zárt kör átvitele
% Zero/pole/gain:
%
% 0.1397 z^2 (z-0.9048) (z-0.9048) (z-0.8187) (z-0.8187) (z-1) (z+0.9048)
% ------------------------------------------------------------------------
% z^2 (z-1) (z-0.9048) (z-0.9048) (z-0.8187) (z-0.8187) (z+0.2656)
%
% (z^2 - 1.266z + 0.4759)
%
%
% Sampling time: 1
step(T)
SMITH-prediktor
két átviteli fv. egyenlőségéből (C hullám és C...):
C~
C = --------------------
1+(1-e^(-sTh))*C~*P+
ez a Smith-prediktor tervezési képlete.
C~
C=--------------------
1+(1-z^(-k))*C~*G+
Példa:
1 1
P+ = ------- * -------
1+10s 1+5s
Gp=... C1=...
L2=C1*Gp % Zero/pole/gain: % 0.0090559 (z-0.9048) (z-0.8187) (z+0.9048) % ------------------------------------------ % z (z-1) (z-0.9048) (z-0.8187) % % Sampling time: 1 [mag,phase,w]=bode(L2); k2=margin(mag,phase-60,w) % k2 = % 31.8911 Ct=k2*C1 % Transfer function: % 31.89 z^2 - 54.96 z + 23.62 % --------------------------- % z^2 - z % Sampling time: 1 margin(Ct*Gp) % a diagrammot egyszercsak levágja. Itt van a határ frekvencia CS= Ct/(1+(1-z^(-1))*Ct*Gp) % Smith-prediktoros szabályzó % Zero/pole/gain: % 31.8911 z^2 (z-0.9048) (z-0.9048) (z-0.8187) (z-0.8187) (z-1) % ------------------------------------------------------------- % z (z-1)^2 (z-0.9048) (z-0.8187) (z^2 + 0.2888z + 0.2613) % % Sampling time: 1 TS=CS*G/(1+CS*G) % Zero/pole/gain: % % 0.2888 z^4 (z-0.9048)^4 (z-0.8187)^4 (z-1)^3 (z+0.9048) (z^2 + 0.2888z + 0.2613) % % --------------------------------------------------------------------------------- % z^3 (z-1)^3 (z-0.9048)^4 (z-0.8187)^4 (z-3.648e-008) (z+3.648e-008) % % (z^2 - 0.7113z + 0.2613) (z^2 + 0.2888z + 0.2613) % % % Sampling time: 1 step(TS,T) % ennek a kettőnek ugyanannak kéne lennie. Valszeg azért nem lett jó, mert egy helyen kézzel írtunk be egy számot 4 tizedesig, ami pontatlan. => majd küld egy mathlab szkriptet
Itt végeztünk az órai anyaggal, nézünk a házihoz egy kis szimulinket.
>> simulink
