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/(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=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
% Táblázatos módszer: 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) - figyelem! _nem_
% dB-ben! - , míg a fázis phase(3) (fokban).
[mag,phase,w]=bode(w0);
% 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. Pl. 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-t 1/mag(mini)-nek választjuk (mag elemei
% _nem dB-ben_ értendők!)
Ap=1/mag(mini);
wc_p=wc_p*Ap % A szabályozó erősítése Ap
w0=wc_p*wp % A felnyitott kör
figure(2); % Rajzolás a 2. grafikus ablakba
margin(w0); % A felnyitott kör Bode-diagramja fázis- és erősítéstartalékkal
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
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
PD szabályzó tervezése
% 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
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
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
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!
% 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 = Ugrásválasz
% 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:
% A vPID(t) zártkörbeli ugrásválaszának maximuma t=0-ban van. Cél hogy ez umax legyen.
% 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];
