Szabályozástechnika - Diszkrétidejű állapotteres szabályozók tervezése

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

 

Simulink modellek

Simulink modellek

1. Töltsd le!
2. Csomagold ki!
3. Másold be a Matlab aktuális munkakönyvtárába!

A mechanikai lengőrendszer leírása

%% Állapotteres szabályozás diszkrét időben

%% Mechanikai lengőrendszer leírása

% A rendszer paraméterei 
m=2;    % A test tömege
k=0.75;  % Rugóállandó
b=0.25; % Csillapítás  

% A folytonos idejű állapotteres leírás mátrixai:

Ac = [0 1; -k/m -b/m]
Bc = [0; 1/m]
Cc = [1 0]
Dc = 0

sys=ss(Ac,Bc,Cc,Dc); % A folytonos idejű rendszer összeállítása

Ts=0.2;                 % Mintavételi idő
dsys=c2d(sys,Ts,'zoh'); % Áttérés diszkrét időre

step(sys,dsys); % A folytonos és diszkrét idejű rendszer ugrásválasza

% A diszkrét idejű állapotteres leírás mátrixai:
[Phi,Gamma,C,D]=ssdata(dsys)
 

Állapotvisszacsatolás tervezése

% Az állapotteres szabályzás alapelve, hogy a zárt körben előre meghatározott pólusokat szeretnénk.
% Ezek általában egy domináns komplex konjugált póluspár (zdom1 és zdom2 - ez határozza meg a tranzienseket),
% és megfelelő számú, a domináns póluspárnál 3-5ször gyorsabb pólus (zcinf).
% Ezekből a segédpólusokból n-2 darabra van szükség, ahol n az állapotváltozók száma.
% Ezt a kívánt póluselrendezést úgy érhetjük el, hogy a bemenetre egy K erősítésvektoron keresztül negatívan
% visszacsatoljuk az állapotváltozókat, azaz u = -K*x. Ha ezt behelyettesítjük a rendszer állapotegyenletébe,
% akkor x[k+1] = Phi*x[k] + Gamma*u[k] = Phi*x[k] + Gamma*(-K*x[k]) = (Phi-Gamma*K)*x[k] lesz a módosult
% rendszer állapotegyenlete. Látható hogy úgy kell a K értékét megválasztani, hogy az Phi-Gamma*K módosult
% rendszermátrix sajátértékei éppen az általunk előírt pólusok legyenek. Az ehhez szükséges K erősítés
% SOR-vektor az Ackermann képlettel egyszerűen meghatározható.
%
% Az állapot-visszacsatolást tartalmazó szabályzó hatásvázlata:

% Gyorsabb, de jobban csillapított zárt kört szeretnénk - az előírásokat
% folytonos időben adjuk meg!
w0=1;
xi=0.8;

% A zárt kör sajátértékei s-ben
% Ha a rendszernek 2-nél több állpotváltozója van, akkor n-2 segédpólust
% (scinf) is fel kell vennünk, melyek 3-5ször gyorsabbak mint a domináns póluspár.
sdom1=-w0*xi+j*w0*sqrt(1-xi^2);
sdom2=conj(sdom1);

% Áttérés z-re az ismert z = exp(s*Ts) képlet alapján
zdom1=exp(sdom1*Ts);
zdom2=exp(sdom2*Ts);
% zcinf=exp(scinf*Ts);

% A zárt kör sajátértékeit tartalmazó vektor, valamint ha n>2,
% akkor n-2 multiplicitással a zcinf pólusok is
phic=[zdom1 zdom2];
% A zárt kör karakterisztikus polinomja
polyphic=poly(phic); 

% Az irányíthatóság ellenőrzése
Mc=ctrb(Phi,Gamma); % Az irányíthatósági mátrix...
rank(Mc)            % ... és rangja
% Ha rank(Mc)=n, akkor a rendszer irányítható

% Állapotvisszacsatolás tervezése az Ackermann-képlet segítségével
K=acker(Phi,Gamma,phic)

% A megfelelő Simulink-modell megnyitása
open('discrete_1');
% Ugyanaz az feladat, mint a folytonos esetben...
% VIGYÁZAT: Minden építőelemnél, ami nem a szakasz része, meg kell
% adni a Ts mintavételi periódusidőt!!!

% Várakozás billentyűlenyomásra
pause
 

Alapjel miatti korrekció

% Amikor az állapot-visszacsatolás 0-ba (alaphelyzetbe) viszi a rendszert, a beavatkozó jel is 0 lesz,
% azaz beáll a stabilis egyensúlyi állapot. Azonban a szabályozásnak nem feltétlenül az a célja,
% hogy 0-ba irányítsunk, hanem célszerű, ha alapjelet is tud követni az eszköz. Ehhez az állapot-visszacsatolót
% "átverjük", az alábbi hatásvázlatnak megfelelően (Nx - oszlopvektor, Nu - skalár):

% Adott a szakaszunk állapotegyenletei: x[k+1] = Phi*x[k] + Gamma*u[k] és y[k] = C*x[k] (legyen D=0)
% Egységugrás alapjel követése esetén célunk, hogy állandósult állapotban a kimenet 1 értékű legyen.
% Továbbá tudjuk, hogy állandósult állapotban, azaz "végtelenben" x[k+1] = x[k]
% A hatásvázlatról látszik, hogy állandósult állapot esetén a K erősítő bemenetén 0 kell, hogy legyen.
% Ekkor x[k] = Nx*r[k] = Nx, mivel r[k]=1 egységugrás alapjel esetén.
% Ha azonban a K bemenete 0, akkor a kimenete is 0, így u[k] = Nu*r[k] = Nu.
% Ezek lapján felírható az alábbi egyenletrendszer:
%
% x[k] = Phi*x[k] + Gamma * u[k]  --> 0 = (Phi-I)*x[k] + u[k]
% 1    = C  *x[k]
%
% Behelyettesítve x[k]=Nx és u[k]=Nu értékeket és az egyenleteket mátrixos alakba átírva:
%
% ( Phi-I   Gamma )   ( Nx )   ( 0 )
% (   C       0   ) * ( Nu ) = ( 1 )
%
% Melyből már kapásból adódik az Nx és Nu erősítések:
%
% ( Nx )   ( Phi-I   Gamma )^-1   ( 0 )
% ( Nu ) = (   C       0   )    * ( 1 )

% A keresett erősítésvektor meghatározása:
% FIGYELEM: Eltérés a folytonos időtől, hogy itt Phi-eye(n) szerepel első elemként.
N=inv([Phi-eye(2) Gamma; C 0])*[0;0;1];
% Az oszlopvektorba n darab nullát kell pakolni és a végére egyetlen 1-est.

% Az Nx-et és Nu-t tartalmazó vektor szétválasztása
Nx=N(1:2) % Annyi elem, ahány állapotunk van (n)
Nu=N(end) % Skalár

% A megfelelő Simulink-modell megnyitása
open('discrete_3');

% Várakozás billentyűlenyomásra
pause
 

Állapotmegfigyelő tervezése

% Az állapot-visszacsatolásnál minden egyes időpillanatban szükségünk van az állapotok aktuális értékeire.
% Ez a gyakorlatban mérésekkel lenne megvalósítható, ám ez nem mindig lehetséges, vagy ha lehetséges,
% akkor csak nagyon drágán. Éppen ezért használunk állapotmegfigyelőt, ami képes minden időpillanatban
% nagyon jó közelítéssel előállítani az állapotok aktuális értékeit, a szakasz kimenetének (y) és
% bemenetének (u) ismeretében.
% A megfigyelő maga is egy lineáris diszkrét idejű rendszer, melynek differenciaegyenlete:
% xhat[k] = F * xhat[k-1] + G * y[k] + H *u[k-1]
% Ahol xhat a becsült állapotok OSZLOP-vektora, valamint dim{xhat}=dim{x}=n.
% Továbbá éljünk az F = Phi-G*C*Phi és H = Gamma-G*C*Gamma választással.
% A becslés hibájára (xtilde = xhat - x) így az alábbi differenciaegyenlet adódik: xtilde[k+1] = F * xtilde[k]
% Látható, hogy az F mátrix sajátértékei fogják meghatározni, hogy milyen gyorsan csengjenek le a megfigyelés
% tranziensei, azaz hogy milyen pontos legyen a megfigyelésünk. Így szeretnénk ha az F mátrix sajátértékei az
% általunk előírt gyorsaságú (a domináns póluspárnál ~5ször gyorsabb) zoinf pólusok lennének.
% phio(z) = det (zI - F) = det (zI-(Phi-G*C*Phi)) = det (zI-(Phi'-Phi'*C'*G')
% A fenti egyenlőség azért igaz, mivel egy rendszer és annak duálisának azonos a karakterisztikus polinomja.
% Így visszavezettük a feladatot egy fiktív rendszerhez történő állapot-visszacsatolás (K2=G') megtervezésére,
% mely fiktív rendszer rendszermátrixa Phi' (Phi transzponált), bemeneti erősítésvektora pedig
% Phi'*C' (Phi transzponált * C transzponált).
% Ezek ismeretében az állapotmegfigyelő G vektora (vagyis annak transzponáltja) már számítható az
% Ackermann képlettel.
%
% Az állapotmegfigyelőt is tartalmazó állapot-visszacsatolásos szabályzó hatásvázlata:

% A megfigyelő sajátértéke s-ben
soinf=-5

% Áttérés z-tartományra
zoinf=exp(soinf*Ts);

% A megfigyelő karakterisztikus polinomjának gyökei - zoinf megfelelő multiplicitással (n)
phio=[zoinf zoinf]
% A megfigyelő karakterisztikus polinomja
polyphio=poly(phio)

% A megfigyelhetőség ellenőrzése
Mo=obsv(Phi,C); % A megfigyelhetőségi mátrix...
rank(Mo)        % ... és rangja
% Ha rank(Mo)=n, akkor a rendszer megfigyelhető!

% Megfigyelő tervezése - VIGYÁZAT: Itt a paraméterezés nem analóg a folytonos esettel!
G=acker(Phi',Phi'*C',phio)'
F=Phi-G*C*Phi
H=Gamma-G*C*Gamma

% A megfelelő Simulink-modell megnyitása
open('discrete_2');
% VIGYÁZAT: Itt nem lehet discrete state space objektumként megadni az állapotmegfigyelőt!
% Elemenként kell azt összeraknunk azt, a differenciaegyenlete alapján!
%
% ^          ^
% x[k] = F * x[k-1] + G * y[k] + H * u[k-1]
%
% FIGYELEM: Mivel aktuális megfigyelőt terveztünk, így az y[k] nincs késleltetve!

% Várakozás billentyűlenyomásra
pause
 

Terhelésbecslő tervezése

% Tegyük fel, hogy a bemenetre rárakódik egy szakaszosan hosszú ideig konstans zaj. Az állapotmegfigyelővel
% ezt is becsülni tudjuk és ezt a becsült értéket kivonhatjuk a bemenetből (beavatkozó jelből), így
% kompenzálhatjuk a zavarás hatását. Ezt úgy oldjuk meg, hogy a zajt egy új állapotváltozónak (xd) tekintjük.
% Mivel a zaj szakaszosan hosszú ideig konstans, ezért xd[k+1]=xd[k], a bemenettől pedig független, ezért a
% differenciaegyenlete: xd[k+1] = 0*x[k] + xd[k] + 0*u[k]. Viszont a többi változó differenciaegyenletébe már
% beleszól az xd zavarást modellező fiktív állapotváltozó, méghozzá a Gamma bemeneti mátrixon keresztül:
% x[k+1] = Phi*x[k] + Gamma*( xd[k]+u[k] )
% Tehát a kibővített rendszerünk állapotegyenletei:
%
% ( x[k+1]  )   ( Phi   Gamma )   ( x[k]  )   ( Gamma )
% ( xd[k+1] ) = (  0      1   ) * ( xd[k] ) + (   0   ) * u[k]
%
%                           ( x[k]  )
%    y[k]     = ( C   0 ) * ( xd[k] )
%
% Új jelöléseket bevezetve a kibővített rendszerünk állapotegyenletei:
%
% xtilde[k+1] = Phitilde * xtilde[k] + Gammatilde * u[k]
%    y[k]     = Ctilde   * xtilde[k]
%
% Tehát az állapotmegfigyelőnket most ehhez a kibővített "tilde" rendszerhez kell megterveznünk.
% A módosult állapotmegfigyelő differenciaegyenlete a hatásvázlatról leolvasható!
% Az terhelésbecslővel kiegészített szabályzó hatásvázlata:

% A kibővített rendszer mátrixai
Phitilde=[Phi Gamma; 0 0 1]; % n darab nulla és fixen 1 darab egyes (SISO)
Gammatilde=[Gamma;0];        % Fixen 1 darab nulla a végére (SISO)
Ctilde=[C 0];                % Fixen 1 darab nulla a végére (SISO)

% A megfigyelő sajátértékeit tartalmazó vektorban zoinf most egyel nagyobb
% multiplicitással szerepel (n+1), hiszen felvettünk egy új (fiktív) állapotot
phiotilde=[zoinf zoinf zoinf];

% Megfigyelőtervezés a kibővített rendszerhez
% Ugyanaz, mint az állapotmegfigyelőnél, csak most a 'tilde rendszerre
Gtilde=acker(Phitilde',Phitilde'*Ctilde',phiotilde)'
Ftilde=Phitilde-Gtilde*Ctilde*Phitilde;
Htilde=Gammatilde-Gtilde*Ctilde*Gammatilde;

% A megfelelő Simulink-modell megnyitása
open('discrete_4');
% t=10 secnél egy egységugrás jellegű zavarás adódik a szakasz bemenetére.
% A modellben a K,Nu és Nx paraméterek ugyanazok, mint amiket korábban meghatároztunk.
% FONTOS: A Simulink alapból 10 secundumig számol, szóval ezt az időt át kell írni 20-ra az ablak tetején!

% Várakozás billentyűlenyomásra
pause
 

Integráló szabályozás tervezése

% Az integráló szabályzó célja a zavarelnyomás és a paraméterbizonytalanságok kiküszöbölése. Ezt úgy érjük el,
% hogy új állapotként felvesszük a kimenet integrálját: xI = integrál (0->t) y(tau) dtau
% Ebből a baloldali téglalap szabályt (LSR) alkalmazva: xI[k+1] = xI[k] + T*y[k] = xI[k] + T*C*x[k]
% Ezzel már felírhatók a kibővített rendszer állapotegyenletei:
%
% ( x[k+1]  )   ( Phi   0 )   ( x[k]  )   ( Gamma ) 
% ( xIfk+1] ) = ( T*C   1 ) * ( xI[k] ) + (   0   ) * u[k]
%
%                             ( x[k]  )
%    y[k]     = ( C   0 )   * ( xI[k] )
%
% Új jelöléseket bevezetve a kibővített rendszerünk állapotegyenletei:
%
% xi[k+1] = Phii * xi[k] + Gammai * u[k]
%  y[k]   = Ci   * xi[k]
%
% Most ehhez a kibővített rendszerhez kell egy új Ktilde = [Kt Ki] állapot-visszacsatolást megterveznünk.
% Az integráló hatást is tartalmazó szabályzó hatásvázlata:

% A kibővített rendszer mátrixai
Phii=[Phi zeros(2,1);C*Ts 1]; % Az első sorban n*1-es nullmátrix
                              % A második sorban fixen 1 darab nulla (SISO)
Gammai=[Gamma;0];             % Fixen 1 darab nulla a végére

% Az integrátor állapotának s = -3 folytonosidejű pólusnak megfelelő sajátértéket írunk elő
phictilde=[zdom1 zdom2 exp(-3*Ts)];

% Állapotvisszacsatolás számítása
Ktilde=acker(Phii,Gammai,phictilde);

% Az állapotvisszacsatolás vektorának felbontása
Kt=Ktilde(1:2); % Annyi eleme van, ahány valódi állapotunk (n)
Ki=Ktilde(3);   % Skalár 

% A megfelelő Simulink-modell megnyitása
open('discrete_5');
% Vigyázat ez itt a terhelésbecslő nélküli modell továbbfejlesztése.
% Az integráló szabályozás is a bemenetre szuperponálódott zavarjelek kiküszöbölésére való.
% Itt Nu helyett egy Ki erősítés van és egy integrátor, valamint K helyett Kt !!!

% Várakozás billentyűlenyomásra
pause