Szabályozástechnika - Soros kompenzátorok tervezése
A VIK Wikiből
A szakasz és a tervezési előírások megadása
%% A szakasz a szokásos:
% Aplant
% Wp(s) = ---------------------
% (1+sT1)(1+sT2)(1+sT3)
Aplant=5;
T1=10;
T2=4;
T3=1;
% A szakasz átviteli függvénye
wp=tf(Aplant,conv(conv([T1 1],[T2 1]),[T3 1]))
% Tervezési előírás: 60 fokos fázistartalék
Phit=60;
P szabályzó tervezése
% Kompenzálási stratégia: Ap segítségével beállítjuk a megfelelő fázistartalékot.
%
% Wc_p(s) = Ap
%
wc_p=tf(1,1) % A szabályzó átviteli függvénye Ap=1 választás mellett
w0=wc_p*wp % A felnyitott kör átviteli függvénye
figure(1); % Rajzolás az 1. grafikus ablakba
bode(w0); % A felnyitott kör Bode-diagramja
% Leolvasás: Azon a körfrekvencián, ahol a felnyitott kör fázisa -120 fok,
% azaz 180+Phit fok, ott a felnyitott kör erősítése 3.62dB. Tehát ahhoz,
% hogy az adott körfrekvencián az erősítés 1 = 0dB legyen (azaz az adott
% frekvencia legyen a vágási körfrekvencia), Ap-t -3.62 dB-nek kell
% választanunk, azaz a körerősítést csökkentenünk kell. Ezzel együtt a
% vágási körfrekvencia csökken, a zárt rendszer is lelassul kissé.
Ap=10^(-3.62/20) % Ap kiszámítása a dB-es értékből
% ---------------------------------------------------------------------------------------------
% Ugyanez számítható egy mechanikus táblázatos módszerrel is:
% A bode függvényt [mag,phase,w]=bode(w0) módon hívva az nem az ábrával, hanem három
% oszlopvektorral tér vissza, amik a frekvencia (w) - amplitúdó (mag) - fázis (phase)
% táblázat oszlopait tartalmazzák. Például a harmadik visszaadott frekvenciaértéken
% (w(3) - a w vektor harmadik eleme) az abszolút érték mag(3), míg a fázis phase(3) (fokban).
% FIGYELEM: A mag oszlopban az magnitúdó értékek NEM dB-ben vannak megadva!!!
[mag,phase,w]=bode(w0);
% [mag,phase,w]=bode(w0,0.01:0.00001:10);
% Így növelhető a pontosság, mivel alapból csak néhány 100 pontban számol a bode függvény.
% Azonban ügyelni kell rá, hogy a keresett frekvencia biztosan benne legyen a megadott
% intervallumban (itt: 0.01-10), azonban túl nagyra sem érdemes választani mert kiakad a Matlab.
% Célszerű előtte megnézni az alap bode(w0) diagramon hogy körülbelül milyen intervallumot kell megadni.
% A teendőnk, hogy megkeressük, melyik visszaadott frekvencián van a fázis
% a legközelebb (akár pozitív, akár negatív irányban) a Phit+180 értékhez.
% A [minv,mini]=min(x) függvény nem csak az x vektor legkisebb elemét
% (minv), hanem annak indexét (mini) is visszaadja.
% Például ha x=[2 1 5], akkor minv=1 és mini=2.
[minv,mini]=min(abs(phase+180-Phit))
% Az erősítés az adott frekvencián mag[mini], mi pedig ott szeretnénk látni
% a vágási körfrekvenciát, így Ap értékét 1/mag(mini)-nek választjuk.
% FIGYELEM: A mag oszlopban az magnitúdó értékek NEM dB-ben vannak megadva!!!
Ap=1/mag(mini);
% ---------------------------------------------------------------------------------------------
wc_p=wc_p*Ap % A szabályozó erősítése Ap szeresére nőtt
w0=wc_p*wp % A felnyitott kör újbóli előállítása a módosult szabályzóval
figure(2); % Rajzolás a 2. grafikus ablakba
margin(w0); % A felnyitott kör Bode-diagramja fázis- és erősítéstartalékkal
% Itt lehet ellenőrizni, hogy sikerült-e elérni a 60 fokos fázistartalékot!
% +/- pár fokos eltérés belefér, de ha nagyobb, akkor valamit elrontottunk!
wcl_p=feedback(w0,tf(1,1),-1) % A zárt kör átviteli függvénye
% Egységnyi merev negatív visszacsatolás esetén tömören:
% wcl_p=feedback(w0,1)
figure(3); % Rajzolás a 3. grafikus ablakba
step(wcl_p); % A zárt kör ugrásválasza
pause; % Várakozás gombnyomásra
close all; % Az összes nyitva lévő grafikus ablak bezárása
PI szabályzó tervezése
% Kompenzálási stratégia:
% 1. Ti-vel kiejtjük a szakasz leglassabb pólusát
% 2. Ap-vel beállítjuk a megfelelő fázistartalékot
%
% ( 1 ) s*Ti + 1
% Wc_pi(s) = Ap * ( 1 + ------ ) = Ap * ----------
% ( s*Ti ) s*Ti
Ti=T1; % A T1-hez (T1 a legnagyobb időállandó) tartozó pólus kiejtése
wc_pi=tf([Ti 1],[Ti 0]) % A szabályzó átviteli függvénye Ap=1 választással
figure(1);
pzmap(wp,'b',wc_pi,'r'); % A pólus-zérus kiejtés szemléltetése
% A fenti hívás wp és wc_pi pólus-zérus elrendezését egy ábrára rajzolja
% fel, a 'b' opció következtében wp pólusait és zérusait kék (blue), míg az
% 'r' opció következtében wc_pi pólusait és zérusait piros (red) színnel.
w0=wc_pi*wp % A felnyitott kör összeállítása
w0=minreal(w0) % A felnyitott kör átviteli függvényének egyszerűsítése
figure(2);
bode(w0); % A felnyitott kör Bode-diagramja
% Leolvasva (természetesen a táblázatos módszer is használható):
Ap=10^(-12.2/20)
wc_pi=Ap*wc_pi % Ap erősítés beállítása a szabályozóban
w0=minreal(wc_pi*wp) % A felnyitott kör (egy lépésben egyszerűsítve is)
figure(3);
margin(w0); % A felnyitott kör Bode-diagramja fázis- és erősítéstartalékkal
wcl_pi=feedback(w0,tf(1,1),-1); % A zárt kör átviteli függvénye
figure(4);
step(wcl_pi); % A zárt kör ugrásválasza
pause; % Várakozás gombnyomásra
close all; % Az összes nyitva lévő grafikus ablak bezárása
Közelítő PD szabályzó tervezése
% Nem ideális deriváló, mert azt nem lehet realizálni, beviszünk egy
% új pólust (-1/Tc), így realizálható lesz.
% Kompenzálási stratégia:
% 1. Beállítjuk az N=Td/Tc arányt
% 2. Td+Tc=(N+1)*Tc-vel kiejtjük a szakasz 2. leglassabb pólusát
% 3. Ap-vel beállítjuk a megfelelő fázistartalékot
%
% ( s*Td ) s*(Td+Tc) + 1 s*(N+1)*Tc + 1
% Wc_pi(s) = Ap * ( 1 + ---------- ) = Ap * --------------- = Ap * ----------------
% ( s*Tc + 1 ) s*Tc + 1 s*Tc + 1
N=10; % Td/Tc arány beállítása - Ökölszabályként N=10
Tc=T2/(N+1); % A T2-höz (2. legnagyobb időállandó) tartozó pólus kiejtése
Td=N*Tc; % Td meghatározása
wc_pd=tf([Td+Tc 1],[Tc 1]) % A szabályozó átviteli függvénye Ap=1 mellett
w0=minreal(wc_pd*wp) % A felnyitott kör átviteli függvénye
figure(1);
pzmap(wp,'b',wc_pd,'r'); % A pólus-zérus kiejtés szemléltetése
figure(2);
bode(w0); % A felnyitott kör Bode-diagramja
% Leolvasva (természetesen a táblázatos módszer is használható):
Ap=10^(2.29/20)
wc_pd=Ap*wc_pd % A szabályozó átviteli függvénye
w0=minreal(wc_pd*wp) % A felnyitott kör átviteli függvénye
figure(3);
margin(w0); % A felnyitott kör Bode-diagramja fázis- és erősítéstartalékkal
wcl_pd=feedback(w0,tf(1,1),-1) % A zárt kör átviteli függvénye
figure(4);
step(wcl_pd); % A zárt kör ugrásválasza
pause; % Várakozás gombnyomásra
close all; % Az összes nyitva lévő grafikus ablak bezárása
Közelítő PID szabályzó tervezése
% Kompenzálási stratégia:
% 1. Beállítjuk az N=Td/Tc arányt
% 2. A szakasz két leglassabb pólusának kiejtése Tau1-el és Tau2-vel, majd
% ezekből a szabályozó Ti, Td és Tc időállandóinak meghatározása
% 3. Ap-vel beállítjuk a megfelelő fázistartalékot
%
% ( 1 s*Td ) Ap s^2*Ti*(Td+Tc) + s*(Ti+Tc) + 1
% Wc_pid(s) = Ap * ( 1 + ------ + ---------- ) = ---- * --------------------------------
% ( s*Ti s*Tc + 1 ) Ti s * (s*Tc + 1)
N=10; % Td/Tc arány meghatározása
Tau1=T1; % Pólus-zérus kiejtés a két leglassabb pólussal
Tau2=T2;
% Az egyenletek felírása Tau1, Tau2 és Ti, Td, Tc kapcsolatára:
% 1. Tau1+Tau2 = Ti+Tc
% 2. Tau1*Tau2 = Ti*(Td+Tc) = Ti*(N+1)*Tc
% Az 1. egyenletből Ti-t kifejezve (Ti=Tau1+Tau2-Tc)
% és T2-be helyettesítve
% egy másodfokú egyenletet kapunk Tc-re, aminek megoldása:
Tcs=roots([N+1 -(N+1)*(Tau1+Tau2) Tau1*Tau2])
% Tc mindenképpen pozitív és valós, ha mindkét gyök az lenne, akkor a
% kisebbet választjuk
Tc=Tcs(2) % Itt a gyökök közül mindkettő valós és pozitív, a 2. a kisebb
% Automatizált megoldás - vájtfülű érdeklődőknek, nem része az anyagnak!
% A find függvény segítségével megkeressük Tcs nulla
% képzetes részű (imag(Tcs)=0) és (&) pozitív valós részű (Tcs>0) elemeit.
% A find függvény a feltételnek eleget tevő elemek indexével tér vissza,
% így az értékek maguk a Tcs(find(...)) módon kaphatók meg. Ezek közül
% pedig a minimálisat választjuk.
Tc=min(Tcs(find(imag(Tcs)==0 & Tcs>0)))
Td=N*Tc % A Td/Tc arány alapján
Ti=Tau1+Tau2-Tc % Az 1. egyenletből adódik
wc_pid=tf(1,1)+tf(1,[Ti 0])+tf([Td 0],[Tc 1]) % A szabályozó Ap=1 mellett
w0=minreal(wc_pid*wp) % A felnyitott kör átviteli függvénye
figure(1);
pzmap(wp,'b',wc_pid,'r'); % A pólus-zérus kiejtés szemléltetése
figure(2);
bode(w0); % A felnyitott kör Bode-diagramja
% Leolvasva (természetesen a táblázatos módszer is használható):
Ap=10^(2.16/20)
wc_pid=Ap*wc_pid % A szabályzó átviteli függvénye a megfelelő Ap mellett
w0=minreal(wc_pid*wp) % A felnyitott kör
figure(3);
margin(w0); % A felnyitott kör Bode-diagramja fázis- és erősítéstartalékkal
wcl_pid=feedback(w0,tf(1,1),-1); % A zárt kör átviteli függvénye
figure(4);
step(wcl_pid); % A zárt kör ugrásválasza
%% A PID szabályzó által kiadott beavatkozó jel
wur_pid=feedback(wc_pid,wp,-1); % Az r->u átviteli függvény zárt körben
% Ilyenkor az előrevezető ág a szabályzó, a visszacsatolás pedig a szakasz!
figure(5);
step(wur_pid); % Az r->u átviteli függvény ugrásválasza
pause; % Várakozás gombnyomásra
close all; % Az összes nyitva lévő grafikus ablak bezárása
A soros kompenzátorok összehasonlítása
figure(1);
% Az összes zárt kör ugrásválasza egy ábrán
step(wcl_p,wcl_pi,wcl_pd,wcl_pid);
legend('P','PI','PD','PID'); % Jelmagyarázat hozzáadása
pause; % Várakozás gombnyomásra
close all; % Az összes nyitva lévő grafikus ablak bezárása
PID szabályzó tervezése maximális beavatkozó jelre
umax=10; % Beavatkozó jel maximális értek (t=0-ban lép fel)
% A szakasz két leglassabb pólusának kiejtése, azonban itt nincs előre
% rögzített N=Td/Tc arány!
% Tudjuk hogy T1+T2=Tc+Ti és hogy T1*T2=Ti*(Td+Tc), ahol T1 és T2 a szakasz
% leglassabb pólusainak időállandói. Ezekből könnyen felírható az alábbi két egyenlet:
% 1. Ti=T1+T2-Tc
% 2. Td=T1*T2/Ti-Tc
% A fenti két egyenletből következik, hogy elég Tc-t paraméterként
% továbbvinnünk, majd ha az megvan, akkor 1. és 2. egyenletekkel adódik Ti és Td.
% Három ismeretlenünk: Ap, Tc, wc (vágási körfrekvencia)
% Három egyenlet:
% I. |W0(j*wc)|-1=0
% II. pi+arg{W0(j*wc)}-Phit=0 FIGYELEM: [Phit]=RADIÁN!!!
% III. vPID(0)-umax=0 vPID = PID szabályzó tag ugrásválasza
% Megoldás az fsolve segítségével, amihez a szükséges függvény
% az UTOLSÓ oldalon került megírásra!!!
% Most jöhet a függvény meghívása!
% De ehhez megfelelő kezdőértékeket kell megadnunk (Tk - 166. oldal)
Ap=umax*(T1+T2-T3)*T3/(T1*T2);
Tc=T3;
wc=1/T3;
x=fsolve('myPID',[Ap,Tc,wc])
Ap=x(1);
Tc=x(2);
% További paraméterek meghatározása a fenti 1. és 2. egyenletekből:
Ti=T1+T2-Tc;
Td=T1*T2/Ti-Tc;
wc_pidu=Ap*(tf(1,1)+tf(1,[Ti 0])+tf([Td 0],[Tc 1]))
w0_pidu=minreal(wc_pidu*wp)
figure(1);
margin(w0_pidu) % Fázistartalék ellenőrzéséhez
wcl_pidu=feedback(w0_pidu,1)
figure(2);
step(wcl_pidu) % Zárt kör ugrásválaszának ellenőrzése
% Maximális beavatkozó jel ellenőrzése
wur_pidu=feedback(wc_pidu,wp)
figure(3);
step(wur_pidu)
A megírandó myPID függvény az fsolve-hoz
% FIGYELEM: Ezt egy külön myPID.m fájlba kell megírni!
% Ez NEM egy általános célú függvény lesz, mivel a T1, T2, Phit és umax
% változókat nem paraméterként kapja meg, hanem az egyszerűség kedvéért
% fixen beleírtam a kódba!!!
% I. egyenlet kifejtése (pólus zérus kiejtést kapásból elvégezve):
% Ap Aplant
% W0(j*wc)= Wc(j*wc)*Wp(j*wc)= --------------------------- * -------------
% Ti * (j*wc) * (1 + j*wc*Tc) (1 + j*wc*T3)
% Ap * Aplant
% |W0(j*wc)|= ------------------------------------------------
% Ti * wc * sqrt{ (1+wc^2*Tc^2) * (1+wc^2*T3^2) }
% Továbbá tudjuk az 1. egyenletből, hogy Ti = T1+T2-Tc,
% valamint a fenti kifejezésből még le kell vonni 1-et hogy megkapjuk f1-et.
% II. Egyenletből f2 már könnyen felírható. A fenti W0(j*wc)-nek tagonként vesszük a fázisát,
% majd ezeket összeadogatjuk (előjelesen). FONTOS, hogy Phit értékét RADIÁNBAN kell megadnunk!!!
% III. Egyenlet kifejtése:
% Amikor a referencia jel megjelenik, a hiba 1, tehát a kezdeti időpillanatban a szabályozó egy 0->1
% átmenetet kap, így a szabályozó válasza megegyezik az egységugrásra adott válasszal (t=0-ban!!).
% Ezután a hiba csökken (ha a szabályozó stabil), ami miatt a beavatkozó jelnek is csökkenni kell.
% Tehát a szabályozó ugrásválaszának t=0 időbeli értéke a beavatkozó jel maximális értékével egyenlő.
% vPID(t=0) = Ap * (1 + Td/Tc)
% Tudjuk a 2. egyenletből, hogy Td=T1*T2/Ti-Tc, ezt beírva a fenti kifejezésbe,
% valamint utána Ti helyére behelyettesítve az 1. egyenletbeli Ti=T1+T2-Tc kifejezést:
% ( T1 * T2 ) T1 * T2
% vPID(t=0) = Ap * ( 1 + ------- - 1 ) = Ap * -----------------
% ( Ti * Tc ) Tc* (T1 +T2 - Tc)
% Ebből levonva umax-ot máris megkaptuk az f3 kifejezést.
function f=myPID(x)
Ap=x(1);
Tc=x(2);
wc=x(3);
Aplant=5;
T1=10;
T2=4;
T3=1;
Phit=60*pi/180;
umax=10;
f1=Ap*Aplant/(T1+T2-Tc)/(wc*sqrt((1+wc^2*Tc^2)*(1+wc^2*T3^2)))-1;
f2=pi-pi/2-atan(wc*T3)-atan(wc*Tc)-Phit;
f3=Ap*T1*T2/(T1+T2-Tc)/Tc-umax;
f=[f1 f2 f3];
