Aktív témák

• #1944 emvy félisten

  2004-11-01 00:38:16 #1944

  Összes hozzászólása itt Válaszok az összes hozzászólására itt Válaszok erre a hozzászólásra

  Privát üzenet küldése

  

  emvy

  félisten

  Bocs srácok, de annyira elegem van, h ezt ki kell engedni itt a fáradt gözben...12 órányi gép előtt ülés...apro kis programocska, a robotos szakirány irányításelmélet házijának első felét képes megoldani, a kunszt az, h univerzálisan, pl. tud megfelelően faktorizálni, stb, szal bárkinek ledokumentálja és megoldja. Kis program, de amíg rájöttem, mi az egésznek az elmélete...szal kemény volt.
  
  disp(' ');
  disp(' ');
  disp(' ');
  disp(' ');
  disp(' ');
  disp('FeedForward-FeedBack szabalyozo-tervezo, identifikalo programocska');
  disp('mate.varga@prohardver.hu ');
  disp(' ');
  % 2 szabadságfokú FF-FB diszkrét szabályozó tervezo algoritmus
  % specifikacios resz
  scr=-3; % sc/abs(re(dominanspoluspar))
  sor=-5; % so/abs(re(dominanspoluspar))
  ksr=0.70; % az eloirt csillapitas erteke
  load meres29; % a betoltendo meresi adatok
  client_ts;
  % inic resz
  s=tf('s');
  z=tf('z',client_ts);
  modeldata=iddata(meres_kimenet__y_,meres_bemenet__u_,client_ts);
  modeldata.InputName='Bemenet';
  modeldata.OutputName='Kimenet';
  
  % identifikacios resz
  ze=modeldata(1:400);
  ze=dtrend(ze);
  domodel=armax(ze,[3 3 2 1]); % diszkret rendszer modellje (DiscreteOpenloopModel)
  comodel=d2c(domodel); % a folytonos rendszer modellje (ClosedOpenloopModel)
  disp(' ');disp('A diszkret rendszer atviteli fuggvenye:')
  domodeltf=tf(domodel.b,domodel.a,client_ts)
  disp(' ');disp('A folytonosideju rendszer atviteli fuggvenye:')
  comodeltf=tf(comodel.b,comodel.a)
  
  % a kért értékek kiszámolása:
  disp('A folytonos rendszer parameterei:');
  comodelb=comodel.b;
  comodelgain=comodelb(4);
  disp(['A:',num2str(comodelgain)]);
  copoles=roots(comodel.a);
  tau=sqrt(1/(copoles(2)*copoles(3)));
  ks=(-1/copoles(2)-1/copoles(3))/(2*tau);
  disp(['tau:',num2str(tau)]);
  disp(['kszi:',num2str(ks)]);
  alfa=-1/copoles(1)/tau;
  disp(['alfa:',num2str(alfa)]);
  disp(' ');
  
  % a referenciamodell ertekeinek szamolasa
  ccpoles=roots([tau*tau 2*ksr*tau 1]);
  ccpoles=roots([tau*tau 2*ksr*tau 1]);
  ccpoles(3)=copoles(1);
  sc=scr*abs(real(ccpoles(1)));
  so=sor*abs(real(ccpoles(1)));
  zc=exp(sc*client_ts);
  zo=exp(so*client_ts);
  drpoles(1)=exp(ccpoles(1)*client_ts); % DiscreteReferencePoles
  drpoles(2)=exp(ccpoles(2)*client_ts);
  drpoles(3)=exp(ccpoles(3)*client_ts);
  disp(['A referenciamodell eloirt polusai diszkret idoben:');
  disp([num2str(drpoles')]);
  disp(' ');
  disp(['A megfigyelo es segedpolusok:');
  disp(['so:',num2str(so),' | sc:',num2str(sc)]);
  disp(['Diszkret idoben ugyanezek:');
  disp(['zo:',num2str(zo),' | zc:',num2str(zc)]);
  
  % faktorizáció, elágazások, miegyéb
  
  % eloszor megallapitjuk a diszkret modell gainjet, az minden esetben kell
  temp=tf(domodel.b,domodel.a,client_ts);
  dozeros=roots(domodel.b);
  domodelzpk=zpk(temp);
  temp=domodelzpk.k;
  dogain=temp(1);
  disp(' ');
  disp('A referenciamodell atviteli fuggvenye:');
  domodelzpk
  disp(' ');
  kiejthetok=0;
  %innentol neha trukkozunk, mert a matlab a specifikaciokkal ellentetben
  %zpkmodellnev.z-re pl. cell-t ad vissza, nem vektort, legalabbis nekem
  %nem sikerult ugy hasznalni, ahogy a doksi mondja.. lasd: 'ltiprops zpk'
  
  % akkor ejtheto ki egy zerus, ha 0 es egy kozotti valos erteke van
  % mivel igy realizalhato folytonos szabalyozoval es stabil szabalyozot
  % eremenyez
  
  % Bminusz (Bm) non monic...
  % tovabba Bp=Bplus
  Bm=dogain;
  Bp=1;
  Bmzeros=[];
  Bpzeros=[];
  
  if ((real(dozeros(2))<1)&(real(dozeros(2))>0))&(isreal(dozeros(2)))
  kiejthetok=1;
  Bp=(z-dozeros(2));
  Bpzeros=dozeros(2);
  else
  Bm=Bm*(z-dozeros(2));
  Bmzeros(1)=dozeros(2);
  end
  if ((real(dozeros(1))<1)&(real(dozeros(1))>0))&(isreal(dozeros(1)))
  kiejthetok=kiejthetok+1;
  Bp=Bp*(z-dozeros(1));
  Bpzeros(2)=dozeros(1);
  
  else
  Bm=Bm*(z-dozeros(1));
  if isempty(Bmzeros)
  Bmzeros=dozeros(1);
  else
  Bmzeros(2)=dozeros(1);
  end
  end;
  
  disp(' ');
  disp('Nem kiejtheto zerusok:');
  disp(num2str(Bmzeros));
  
  
  disp('Kiejtheto zerusok:');
  disp(num2str(Bpzeros));
  
  % polinomfokszamok megallapitasa
  
  bpzsize=size(Bpzeros);
  bpzsize=bpzsize(2);
  bmzsize=2-bpzsize;
  
  disp(' ');
  disp('Polinomok osszeallitasa:');
  disp(' ');
  disp('1 integratort feltetelezunk');
  disp(' ');
  disp('Ao polinom:');
  Ao=(z-zo)^3;
  Ao
  disp(' ');
  disp('Am polinom fokszamfeltetele: Am=2, ha gr B-=0, Am=1+gr B- egyebkent');
  disp(' ');
  if bmzsize==2
  disp('A polinom harmadfoku, tehat zc is polus lesz.');
  Am=(z-drpoles(1))*(z-drpoles(2))*(z-zc);
  Amzeros=[drpoles(1);drpoles(2);zc];
  disp(['Am:');
  Am
  disp(' ');
  disp(['Polusok:');
  disp(num2str(Amzeros));
  else
  disp('A polinom masodfoku, csak az eredeti dominans poluspar szerepel.');
  Am=(z-drpoles(1))*(z-drpoles(2));
  Amzeros=[drpoles(1);drpoles(2)];
  disp(['Am:');
  Am
  disp(' ');
  disp(['Polusok:');
  disp(num2str(Amzeros));
  end
  disp(' ');
  Amsize=size(Amzeros);
  Amsize=Amsize(1);
  Aosize=3;
  Sgr=3;
  disp('gr S=gr A=3');
  disp('gr Ao=gr A=3');
  disp(['gr R1''=gr B- =',num2str(bmzsize)]);
  Rgr=bmzsize;
  disp(' ');
  disp('C=Am*Ao=');
  Am*Ao
  disp(' ');
  [Cnum,Cden]=zp2tf([Amzeros;zo;zo;zo],0,1);
  disp(' ');
  disp('A:=A(z-1)=');
  [Anum,Aden]=zp2tf([drpoles';1],0,1);
  A=tf(Anum,1,1)
  disp(' ');
  disp('B:=B- =');
  [Bnum,Bden]=zp2tf(Bmzeros',0,dogain);
  B=tf(Bnum,1,1)
  disp(' ');
  
  % A diophantoszi egyenletrendszer megoldasanak matrixai
  
  disp('Osszeallitom a linearis egyenletrendszer matrixait');
  if bmzsize==1
  AA=[[Anum',[Bnum 0 0 0]',[0 Bnum 0 0]',[0 0 Bnum 0]',[0 0 0 Bnum]';
  CC=(Cnum(2:6)-[Anum(2:5),0])';
  else
  AA=[[Anum';0],[0;Anum',[Bnum 0 0 0]',[0 Bnum 0 0]',[0 0 Bnum 0]',[0 0 0 Bnum]';
  CC=(Cnum(2:7)-[Anum(2:5),0,0])';
  end
  disp(' ');
  disp('A mátrix:');
  disp(num2str(AA));
  disp(' ');
  disp('C mátrix:');
  disp(num2str(CC));
  disp(' ');
  disp('Megoldva az egyenletrendszert, a paramétervektor:');
  XX=inv(AA)*CC;
  disp(num2str(XX));
  Rparam=XX(1:bmzsize);
  Sparam=XX((bmzsize+1)Amsize+3));
  
  % Celegyenes...
  %Rz' es a tobbi... definicio szerint
  
  Rzv=tf([1;Rparam]',1,client_ts); % 1 kell elore, mert Rzv monic
  Sz=tf(Sparam',1,client_ts); % ez meg nem monic
  Amnum=zp2tf(Amzeros,0,1); % Bm' meghatarozasahoz kellenek az alabbiak
  Bmnum=zp2tf(Bmzeros',0,dogain);
  Bmv=sum(Amnum)/sum(Bmnum);
  disp(' ');
  disp('R1'' polinom:');
  Rzv
  disp(' ');
  Tz=Bmv*Ao;
  Rz=Bp*(z-1)*Rzv;
  FF=Tz/Rz;
  FB=Sz/Rz;
  disp(' ');
  disp('Vegul a szabalyozoban szereplo polinomok:');
  disp(' ');
  disp('T:');
  Tz
  disp(' ');
  disp('S:');
  Sz
  disp(' ');
  disp('R:')
  Rz
  disp(' ');
  disp('Az elorecsatolo ag atviteli fuggvenye az FF, a visszacsatolo ag');
  disp('atviteli fuggvenye az FB valtozoban talalhatok az esetleges tovabbi felhasznalashoz.');
  disp(' ');
  
  
  % Szimulaciok
  disp('A szimulacio a diszkretideju modellt hasznalja. Ellenorzeskepp');
  disp('Simulink hasznalhato a folytonosideju modell szabalyozasahoz,');
  disp('de korrekt identifikacio utan a kimenetek meg fognak egyezni.');
  disp(' ');
  % normal idoallando
  
  sys=feedback(domodelzpk,FB);
  sys=FF*sys;
  sys2=sys;
  y100=step(sys,15);
  sys=feedback(1,domodelzpk*FB); % a vezerlojel megallapitasahoz csak a vissza
  %csatolas modjan valtoztatunk termeszetesen
  sys=FF*sys;
  u100=step(sys,15);
  
  comodelzpk=zpk(comodel);
  comodelgain=comodelzpk.k(1);
  cozeros=roots(comodel.b);
  copoles=roots(comodel.a);
  
  [tempn,tempd]=zp2tf(1,copoles'/1.25,1);
  comodel125=tf(comodelb(4)/1.25/1.25,tempd); % ketszer leosztunk, az elozo sor miatt
  [tempn,tempd]=zp2tf(1,copoles'/0.75,1);
  comodel75=tf(comodelb(4)/0.75/0.75,tempd);
  
  domodelzpk=c2d(comodel125,client_ts);
  sys=feedback(domodelzpk,FB);
  sys=FF*sys;
  y125=step(sys,15);
  sys=feedback(1,domodelzpk*FB); % a vezerlojel megallapitasahoz csak a vissza
  %csatolas modjan valtoztatunk termeszetesen
  sys=FF*sys;
  u125=step(sys,15);
  
  domodelzpk=c2d(comodel75,client_ts);
  sys=feedback(domodelzpk,FB);
  sys=FF*sys;
  y75=step(sys,15);
  sys=feedback(1,domodelzpk*FB); % a vezerlojel megallapitasahoz csak a vissza
  %csatolas modjan valtoztatunk termeszetesen
  sys=FF*sys;
  u75=step(sys,15);
  
  pltime=(1:150)';
  disp(' ');
  disp('Idoallandok: kek: eredeti, zold:125%-os, piros: 75%-os');
  figure(1);
  plot(pltime,y100,pltime,y125,pltime,y75);
  grid on
  figure(2);
  plot(pltime,u100,pltime,u125,pltime,u75);
  grid on
  
  % es kesz is.

Aktív témák