Szabályozástechnika - Soros kompenzátorok tervezése

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

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_pd(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)
%
%
%                    Aplant                                  Aplant 
% Wp(s) = ---------------------------- =  ----------------------------------------------
%         (1+s*Tau1)(1+s*Tau2)(1+s*T3)    (s^2*Tau1*Tau2 + s*(Tau1+Tau2) + 1) * (s*T3+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,
% aztán egy kicsit rendezgetve 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 kell hogy legyen! Ha mindkét gyökre teljesül ez a
% feltétel, akkor a kisebbiket választjuk mindig.

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 p-z kiejtéssel
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

% A P szabályzó relatív gyors, de nagy a túllövése és nagy a maradó hiba.
% A PI szabályzó lassabb, de kicsi a túllövése és zérus a maradó hiba.
% A közelítő PD nagyon gyors, de relatív nagy a túllövése és a maradó hibája.
% A közelítő PID miden egyesíti az előzőek jó tulajdonságait: Gyors, kis túllövés és zérus maradó hiba. 

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 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)

pause; % Várakozás gombnyomásra
close all; % Az összes nyitva lévő grafikus ablak bezárása
 

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];