Szabályozástechnika - 2DOF szabályzó tervezése

A szakasz megadása és a zárt körre vonatkozó előírások

%% A szakasz a szokásos: Aplant / (1+sT1)(1+sT2)(1+sT3)

Aplant=5;
T1=10;
T2=4;
T3=1;

wp=tf(Aplant,conv(conv([T1 1],[T2 1]),[T3 1]))

%% A zárt körre vonatkozó előírások

xi=0.7; % A domináns póluspár csillapítása
w0=1/3; % A domináns póluspár sajátfrekvenciája

scinf=-1; % További pólusok
soinf=-5; % Megfigyelő pólusok

lint=1; % Zárt körbe veendő integrátorok száma

%% A zárt kör domináns póluspárjának számítása

sdom1=-w0*xi+j*w0*sqrt(1-xi^2)
sdom2=conj(sdom1)

%% Áttérés diszkrét időre

Ts=0.2; % Mintavételi idő [s]
dp=c2d(wp,Ts,'zoh') % Áttérés diszkrét időre nulladrendű tartószervet feltételezve
zdom1=exp(sdom1*Ts) % Áttérés a z=exp(s*Ts) képlet szerint
zdom2=exp(sdom2*Ts)
zcinf=exp(scinf*Ts)
zoinf=exp(soinf*Ts)
 

A számláló faktorizálása

B=dp.num{1}  % Számláló kiolvasása a tf struktúrából
A=dp.den{1}; % Nevező kiolvasása a tf struktúrából

% Mivel a polinomok fokszámait majd az őket tároló vektorok 
% hosszából számítjuk (fokszám = vektor elemszáma - 1), a Matlab pedig
% az átviteli függvény számlálóját és nevezőjét egyforma elemszámú 
% vektorokban tárolja, ezért B elejéről le kell választanunk a vezető
% 0 elemeket. Jelen esetben csak 1 darab van, így csak azt kell levágni.
B=B(2:end)
% Általános megoldás - csak ínyenceknek!
B=B(B~=0); % Bminus-ba az eredeti Bminus nem 0 elemeit töltjük

Bpoles=roots(B) % B gyökeinek számítása
% Kiejthető gyökök: Egységkörön belül és NEM tisztán negatív valósak
% Nem kiejthető gyökök: Egységkörön kívűl, vagy tisztán negatív valósak
Bminus=B % Esetünkben egyik zérus sem kiejthető, mivel mind negatív valósak
Bplus=1

% Általánosan kis manuális beavatkozással:
% roots(B)
% Megnézzük hogy mik a kiejthető és mik a ki nem ejthető gyökök
% Bplus=poly([ kiejthető gyökök felsorolása ])
% Bminus=B(1)*poly([ nem kiejthető gyökök felsorolása ])

 
%% Csak ÍNYENCEKNEK - általános megoldás a Matlab lehetőségeit kihasználva.
% Bminus-ba azok a gyökök kerülnek, amiknek az abszolútértéke >=1 ill.
% tisztán valósak és negatívak. Az & jel a két feltétel közötti ÉS
% kapcsolatot jelenti, az union két halmaz (esetünkben vektor) uniója. Ezt a
% műveletet használva minden olyan elem bekerül Bminuspoles-ba, amelyre
% valamelyik feltétel igaz, de akkor is csak egyszer fog szerepelni, ha
% mindkét feltétel igaz rá.
Bminuspoles=union(Bpoles(abs(Bpoles)>=1),...
Bpoles((real(Bpoles)<0)&(imag(Bpoles)==0))); 
% Bplus-ba Bpoles es Bminuspoles különbsége kerül
Bpluspoles=setdiff(Bpoles,Bminuspoles); 
% Bplus polinom összeállítása a gyökökből
Bplus=poly(Bpluspoles); 
% Hogy Bplus*Bminus=B legyen, Bminus-t még szorozni kell egy megfelelő 
% konstanssal, hiszen a poly mindenképpen 1 vezető együtthatójú polinomot
% képez az elemekből és az elmélet szerint Bplus monic (1 vezető együtthatójú).
Bminus=poly(Bminuspoles)*B(1);
 

Fokszámfeltételek (Tk - 240. oldal)

% A fokszámok természetesen kézzel is számíthatók, de itt kihasználjuk a
% Matlab lehetőségeit. 
grBminus=length(Bminus)-1 % A vezető 0-t már korábban leválasztottuk
grA=length(A)-1

% grAm=1+grBminus+{1: grBminus=0; 1 egyébként} - grAm legalabb másodfokú
% kell legyen! 
% grAo=grA+lint-1-{1: grBminus=0; 1 egyébként}
% Itt grBminus=2, így nem kell hozzáadnunk illetve levonnunk semmit.
grAm=1+grBminus
grAo=grA+lint-1

% Általános megoldás - csak ínyenceknek!
% grBminus==0 1-re értékelődik ki ha Bminus fokszáma 0, 0-ra egyébként
grAm=1+grBminus+(grBminus==0); 
grS=grA+lint-1-(grBminus==0)

grR1prime=grBminus

%% Referencia- és megfigyelő polinom számítása
% A polinomokat a gyökeikből a poly függvény segítségével számítjuk.
% Fontos, hogy a poly 1 vezető együtthatójú (monic) polinomokat képez, és
% az elmélet szerint Am és Ao is monic.
% Am harmadfokú (grAm=3), így a domináns póluspár mellett grAm-2, azaz 
% egyszeres multiplicitással kell a zcinf pólust szerepeltetni
Am=poly([zdom1 zdom2 zcinf])
% Általános megoldás
Am=poly([zdom1 zdom2 ones(1,grAm-2)*zcinf]) 

% Ao egyetlen gyöke zoinf, ezt kell most grAo, azaz háromszoros 
% multiplicitással szerepeltetni
Ao=poly([zoinf zoinf zoinf])
% Általános megoldás
Ao=poly(ones(1,grAo)*zoinf) 


%% B'm szamitasa
% B'm = Am(z=1) / Bminus(z=1)
Bmprime=polyval(Am,1)/polyval(Bminus,1)
 

A diophantoszi polinomegyenlet megoldása

polyint=[1 -1]; % Az integrátor polinomja (z-1)

% Általános megoldás - csak ínyenceknek!
% lint számú integrátorra a polyint polinom számítása ( (z-1)^lint )
polyint=1;
for ii=1:lint
    polyint=conv(polyint,[1 -1]);
end;

polyA=conv(A,polyint) % A diophantoszi egyenletben szereplő A polinom
polyB=Bminus % A diophantoszi egyenletben szereplő B polinom
polyC=conv(Am,Ao) % A diophantoszi egyenletben szereplő C polinom

% A polinomegyenletnek megfelelő dioAx=dioB lineáris egyenletrendszer mátrixainak
% összeállítása. A toeplitz(C,R) függvény olyan Toeplitz-mátrix-szal tér
% vissza, melynek első oszlopa C, első sora pedig R. Ha C(1) és R(1) 
% nem egyezik meg, akkor figyelmeztetés mellett C(1) értéke lesz a mátrix
% (1,1) indexű eleme.

% Esetünkben az dioA mátrix két blokkból áll. Mivel az ismeretlenek száma 6,
% (R és S polinomok ismeretlenjei. Polinom ismeretlenjeinek száma = fokszám +1.
% Mivel R monic polinom (legmagyasabb fokszámú tagjának a szorzója 1-es),
% így hiába másodfokú, csak 2 ismeretlen van benne.)
% Ezért mindkét blokk sorainak száma 6, az első blokk oszlopainak száma 2
% (S ismeretlenjeinek száma), a másodiké pedig 4 (R ismeretlenjeinek a száma). 
% A toeplitz blokkok paramétereinek megadásaikor, így megfelelő
% számú 0-t kell beszúrnunk az oszlopban (polyA illetve Bminus) szereplő 
% polinomok után úgy, hogy azok elemszáma 6 legyen. Az első sorok első
% eleme Am(1) illetve Bminus(1), ezek után még 1 illetve 3 darab 0-t kell
% szerepeltetnünk.

% Hogyan is néz ez ki?

% (1	  0	| b0	  0	  0	  0 )   (r1)   (c1-a1)
% (a1	  1	| b1	  b0	  0	  0 )   (r2)   (c2-a2)
% (a2	  a1	| b2	  b1	  b0	  0 )   (s0)   (c3-a3)	
% (a3	  a2	| 0	  b2	  b1      b0) * (s1) = (c4-a4)
% (a4	  a3	| 0	  0	  b2      b1)   (s2)   (c5   )
% (0	  a4	| 0	  0	  0	  b2)   (s3)   (c6   )

dioA=[toeplitz([polyA 0],[polyA(1) 0])...
toeplitz([Bminus 0 0 0],[Bminus(1) 0 0 0])]

% Általános megoldás - csak ínyenceknek!
dioA=[toeplitz([polyA zeros(1,length(polyC)-length(polyA)-1)],...
[1 zeros(1,grR1prime-1)])...
toeplitz([polyB zeros(1,length(polyC)-length(polyB)-1)],...
[polyB(1) zeros(1,grS)])];

% Most a jobb oldalon álló vektor első négy eleme polyC(2)-polyA(2) ...
% polyC(5)-polyA(5) (polyA és polyC első eleme 1, amit nem használunk fel),
% a további sorokban pedig polyC(6) és polyC(7) áll, mivel polyA egy 
% negyedfokú polinom, így nincsenek további együtthatói. A legegyszerűbb 
% ezt a vektort olyan módon számítani, hogy polyC 2..7-ik eleméből kivonunk
% egy olyan vektort, ami polyA 2..5-ik eleme után két 0-t tartalmaz. 
% Ügyelni kell arra, hogy a polinomokat sorvektorként ábrázoljuk, azonban
% a lineáris egyenletrendszer jobb oldalán oszlopvektor áll, így az
% eredményt transzponálnunk kell.

dioB=[polyC(2:end)-[polyA(2:end) 0 0]]'

 
% Általános megoldás - csak ínyenceknek!
% A polyA elemei után beszúrandó 0-k számát a két polinom fokszámkülönbsége
% (a vektorok hosszának különbsége) alapján határozzuk meg.
dioB=polyC(2:end)'-[polyA(2:end) zeros(1,length(polyC)-length(polyA))]';

% Az egyenletrendszer megoldása:
dioSol=inv(dioA)*dioB
 

A szabályzó átviteli függvényék polinomjainak meghatározása

% A dioSol vektor most a következő elemeket tartalmazza:
% dioSol=[r1 r2 s0 s1 s2 s3]'
% Ne feledjük, hogy a megoldás is oszlopvektor formában adódott!

% R1' számítása
% R1' másodfokú, ám mivel monic, ezért legnagyobb fokszámú együtthatója 1,
% így csak a további 2 együttható szerepel a megoldásvektorban. R1prime
% vektor polinomnak felel meg, így sorvektornak kell lennie, úgyhogy
% transzponálnunk kell a megoldásvektor elemeit
R1prime=[1 dioSol(1:2)'] 
% Általános megoldás - csak ínyenceknek!
R1prime=[1 dioSol(1:grR1prime)'] ;

% R polinom számítása
% R=Bplus*(z-1)^l*R1'
R=conv(conv(Bplus,polyint),R1prime)

% S polinom számítása
% S együtthatóit a megoldásvektor 3..6-ik elemei tartalmazzák, de a vektort
% még transzponálnunk kell, hogy S vektor sorvektor legyen
S=dioSol(3:6)'
% Általános megoldás - csak ínyenceknek!
S=dioSol(grR1prime+1:end)';

% T polinom számítása
T=Bmprime*Ao % Bmprime egy skalár, így nincs szükség a conv-ra
 

A zárt szabályozási kör összeállítása

%% A zárt szabályozási kör összeállítása

d_ff=tf(T,R,Ts) % Szabályozó az előrevezető ágban
d_fb=tf(S,R,Ts) % Szabályozó a visszacsatoló ágban
d_cl=d_ff*feedback(dp,d_fb,-1); % A zárt kör számítása
d_cl=minreal(d_cl) % PZ-kiejtés

figure(1);
step(d_cl); % A zárt kör ugrásválasza
xlabel('t');
ylabel('y');

figure(2);
pzmap(d_cl); % A csillapítás és frekvencia leolvasható és ellenőrizhető

% A beavatkozó jelre vonatkozó átviteli függvény zárt körben
d_u=d_ff*feedback(tf(1,1),dp*d_fb); 
figure(3);
step(d_u); % A beavatkozó jel egységugrás alapjel esetén
xlabel('t');
ylabel('u');
 

A robosztusság illusztrációja

% A szakasz minden paraméterét 25%-al megnöveljük ill. lecsökkentjük
wp125=tf(Aplant,conv(conv([1.25*T1 1],[1.25*T2 1]),[1.25*T3 1]));
wp75=tf(Aplant,conv(conv([0.75*T1 1],[0.75*T2 1]),[0.75*T3 1]));
dp125=c2d(wp125,Ts);
dp75=c2d(wp75,Ts);

% Zárt körök összeállítása a perturbalt szakaszokból es az _eredeti_
% szabályozóból
d_cl75=series(d_ff,feedback(dp75,d_fb,-1));
d_cl125=series(d_ff,feedback(dp125,d_fb,-1));

% A beavatkozó jelre vonatkozó átviteli függvények 
d_u75=series(d_ff,feedback(tf(1,1),dp75*d_fb));
d_u125=series(d_ff,feedback(tf(1,1),dp125*d_fb));

% Az ugrásválaszok ábrázolása
figure(4);
step(d_cl,d_cl75,d_cl125);
legend('Nominalis rendszer','75%','125%');
xlabel('t');
ylabel('y');

figure(5);
step(d_u,d_u75,d_u125);
legend('Nominalis rendszer','75%','125%');
xlabel('t');
ylabel('u');