function [reg1,fil1,reg2,fil2,reg3,fil3]=caspa
%CASPA    Menuegefuehrte Koeffiz.-ber. von max. drei kaskadierten dig. Regelern
%         [reg1,fil1,reg2,fil2,reg3,fil3]=caspa
% vom PID-Typ nach dem Digitalen Betragsoptimum; entspricht Mehrfachanwendung
% von bod.m; Programmsteuerung menuegefuehrt
% "caspa" errechnet
%    reg1.. Reglerkoeffizienten des inneren Kaskadenreglers
%    fil1.. Filterkoeffizienten des inneren Kaskadenreglers
%    reg2.. Reglerkoeffizienten des mittleren Kaskadenreglers
%    fil2.. Filterkoeffizienten des mittleren Kaskadenreglers
%    reg3.. Reglerkoeffizienten des äußeren Kaskadenreglers
%    fil3.. Filterkoeffizienten des äußeren Kaskadenreglers
% nach menügesteuerter Parametereingabe der
% Regelstrecke (3.Schleife optional) bzw. der Reglerberechnung
% und errechnet wahlweise Regelabweichung, Über- und Unter-
% schwingweite und eine diskrete Sprungantwort

% SCHLUESSELWORTE:  Digitales Betragsoptimum, Kaskadenregelstruktur,
%                   Programmsteuerung menuegefuehrt, Regleroptimierung,
%                   Digitaler Regler von PID-Typ

% Vorraussetzung ist das Vorhandensein
% der m-Files habod.m und trans.m

% Copyright (c) 1994 Hagen Jahn
%                    TU Dresden; Elektrotechnisches Institut
%                    Mommsenstr. 13; D-01062 Dresden; Germany

% Streckenparametereigabe
% Innere Schleife
disp(' ');
disp(' Kaskadenregelstruktur mit digitaler Regelung');
input('Zwei- (bel.) oder dreischleifigen (3) Struktur ?  ');
if ans==3 sch=3; else sch=0; end
disp(' ');
disp(' Kopplung zwischen Regler und Strecke: ');
disp(' Halteglied-0.Ordng (bel.)/Abtaster (1) =');
k=input('k= '); if k==1  k=0; else k=1; end
disp(' ');
disp(' Geben Sie die Abtastzeit in Milisekunden ein:  ');
T=input('T = ');
disp(' ');
disp(' Folgende Abkuerzungen wurden gewaehlt:');
disp(' ');
disp(' [1] 1 V-A ...1 Verzoegerungsglied mit Augenblickswertmessung');
disp(' [2] 1 V-M ...1 Verzoegerungsglied mit Mittelwertmessung');
disp(' [3] 2 V-A ...2 Verzoegerungsglieder mit Augenblickswertmessung');
disp(' [4] 2 V-M ...2 Verzoegerungsglieder mit Mittelwertmessung');
disp(' [5] 3 V-A ...3 Verzoegerungsglieder mit Augenblickswertmessung');
disp(' [6] 3 V-M ...3 Verzoegerungsglieder mit Mittelwertmessung');
disp(' [7] 1 I-A ...1 Integralgliedglied mit Augenblickswertmessung');
disp(' [8] 1 I-M ...1 Integralgliedglied mit Mittelwertmessung');
disp(' [9]  IV-A ...1 Verzoegerungs- u.1 I-Glied mit Augenblickswertmessung');
disp(' [10] IV-M ...1 Verzoegerungs- u.1 I-Glied mit Mittelwertmessung');
disp(' [11] I2V-A ..2 Verzoegerungs- u.1 I-Glied mit Augenblickswertmessung');
disp(' [12] I2V-M ..2 Verzoegerungs- u.1 I-Glied mit Mittelwertmessung');
ss=input('Waehlen Sie die innere Streckenstruktur aus : ');
disp(' '); co1=0;
if ((ss<1)|(ss>12)|(rem(ss,1)~=0)) error('Falsche Eingabe'); end
if ((ss==1)|(ss==3)|(ss==5)|(ss==7)|(ss==9)|(ss==11))&k==0 
     code1=0;
elseif ((ss==1)|(ss==3)|(ss==5)|(ss==7)|(ss==9)|(ss==11))&k==1 
     code1=1;
elseif ((ss==2)|(ss==4)|(ss==6)|(ss==8)|(ss==10)|(ss==12))&k==0 
     code1=1; co1=3;
else code1=2; end
disp(' Beachtung der Rechen-, Mess- oder Stellgliedreaktionszeit?') 
input('nein(bel.)/ja(1):');
if ans==1
  disp(' '); 
  disp(' Geben Sie die gewuenschte Zeitverschiebung in ms ein!');
  m1=input('m1 = ');
else m1=0; end
disp(' ');
if (ss==1)|(ss==2)
  disp(' Geben Sie die Konstanten nach der Form Vs/(1+pT1) ein ');
  n01=input('Vs = '); T1=input('T1 = '); d01=[T1 1];
elseif (ss==3)|(ss==4)
  disp(' Geben Sie die Konstanten nach der Form Vs/[(1+pT1)(1+pT2)] ein ');
  n01=input('Vs = '); T1=input('T1 = '); T2=input('T2 = ');
  d01=[(T1*T2) (T1+T2) 1];
elseif (ss==5)|(ss==6)	
  disp(' Geben Sie die Konstanten nach der Form Vs/[(1+pT1)(1+pT2)(1+pT3] ein ');
  n01=input('Vs = '); T1=input('T1 = '); T2=input('T2 = '); T3=input('T3 = ');
  d01=[T1*T2*T3 T3*(T1+T2)+(T1*T2) T3+T2+T1 1];
elseif (ss==7)|(ss==8)
  disp(' Geben Sie die Konstanten nach der Form Vs/(pT1) ein ');
  n01=input('Vs = '); T1=input('T1= '); d01=[T1 0];
elseif (ss==9)|(ss==10)
  disp(' Geben Sie die Konstanten nach der Form Vs/[p(p+T1)] ein ');
  n01=input('Vs = '); T1=input('T1= '); d01=[1 T1 0];
else
  disp(' Geben Sie die Konstanten nach der Form Vs/[p(p+T1)(p+T2)] ein ');
  n01=input('Vs = '); T1=input('T1= '); T2=input('T2= ');
  d01=[1 (T1+T2) (T1*T2) 0];
end;
disp(' ');
% Mittlere Schleife
disp(' [1] fuer 1 V-A,[2] fuer 1 V-M,[3] fuer 2 V-A,[4] fuer 2 V-M,');
disp(' [5] fuer 3 V-A,[6] fuer 3 V-M,[7] fuer 1 I-A,[8] fuer 1 I-M');
disp(' [9] fuer IV-A,[10] fuer IV-M,[11] fuer I2V-A,[12] fuer I2V-M');
ss=input('Waehlen Sie die mittlere Streckenstruktur aus : ');
disp(' '); co2=0;
if ((ss<1)|(ss>12)|(rem(ss,1)~=0)) error('Falsche Eingabe'); end
if ((ss==1)|(ss==3)|(ss==5)|(ss==7)|(ss==9)|(ss==11))&k==0
      code2=0;
elseif ((ss==1)|(ss==3)|(ss==5)|(ss==7)|(ss==9)|(ss==11))&k==1
      code2=1;
elseif ((ss==2)|(ss==4)|(ss==6)|(ss==8)|(ss==10)|(ss==12))&k==0
      code2=1; co2=3;
else  code2=2; end
disp(' Beachtung der Rechen-, Mess- oder Stellgliedreaktionszeit?') 
input('nein(bel.)/ja(1):');
if ans==1
  disp(' ');
  disp(' Geben Sie die gewuenschte Zeitverschiebung in ms ein!');
  m2=input('m2 = ');
else m2=0; end
disp(' ');
if (ss==1)|(ss==2)
  disp(' Geben Sie die Konstanten nach der Form Vs/(1+pT1) ein ');
  n02=input('Vs = '); T1=input('T1 = '); d02=[T1 1];
elseif (ss==3)|(ss==4)
  disp(' Geben Sie die Konstanten nach der Form Vs/[(1+pT1)(1+pT2)] ein ');
  n02=input('Vs = '); T1=input('T1 = '); T2=input('T2 = ');
  d02=[(T1*T2) (T1+T2) 1];
elseif (ss==5)|(ss==6)	
  disp(' Geben Sie die Konstanten nach der Form Vs/[(1+pT1)(1+pT2)(1+pT3] ein ');
  n02=input('Vs = '); T1=input('T1 = '); T2=input('T2 = '); T3=input('T3 = ');
  d02=[T1*T2*T3 T3*(T1+T2)+(T1*T2) T3+T2+T1 1];
elseif (ss==7)|(ss==8)
  disp(' Geben Sie die Konstanten nach der Form Vs/(pT1) ein ');
  n02=input('Vs = '); T1=input('T1= '); d02=[T1 0];
elseif (ss==9)|(ss==10)
  disp(' Geben Sie die Konstanten nach der Form Vs/[p(p+T1)] ein ');
  n02=input('Vs = '); T1=input('T1= '); d02=[1 T1 0];
else
  disp(' Geben Sie die Konstanten nach der Form Vs/[p(p+T1)(p+T2)] ein ');
  n02=input('Vs = '); T1=input('T1= '); T2=input('T2= ');
  d02=[1 (T1+T2) (T1*T2) 0];
end;
disp(' ');
if sch==3
  % Aeussere Schleife
  disp(' [1] fuer 1 V-A,[2] fuer 1 V-M,[3] fuer 1 V-A,[4] fuer 2 V-M');
  disp(' [5] fuer 3 V-A,[6] fuer 3 V-M,[7] fuer 1 I-A,[8] fuer 1 I-M');
  disp(' [9] fuer IV-A,[10] fuer IV-M,[11] fuer I2V-A,[12] fuer I2V-M');
  ss=input('Waehlen Sie die aeussere Streckenstruktur aus : ');
  disp(' '); co3=0;
  if ((ss<1)|(ss>12)|(rem(ss,1)~=0)) error('Falsche Eingabe'); end
  if ((ss==1)|(ss==3)|(ss==5)|(ss==7)|(ss==9)|(ss==11))&k==0
        code3=0;
  elseif ((ss==1)|(ss==3)|(ss==5)|(ss==7)|(ss==9)|(ss==11))&k==1
        code3=1;
  elseif ((ss==2)|(ss==4)|(ss==6)|(ss==8)|(ss==10)|(ss==12))&k==0
        code3=1; co3=3;
  else	code3=2; end
  disp(' Beachtung der Rechen-, Mess- oder Stellgliedreaktionszeit?') 
  input('nein(bel.)/ja(1):');
  if ans==1
    disp(' ');
    disp(' Geben Sie die gewuenschte Zeitverschiebung in ms ein!');
    m3=input('m3 = ');
  else m3=0; end
  disp(' ');
  if (ss==1)|(ss==2)
    disp(' Geben Sie die Konstanten nach der Form Vs/(1+pT1) ein ');
    n03=input('Vs = '); T1=input('T1 = '); d03=[T1 1];
  elseif (ss==3)|(ss==4)
    disp(' Geben Sie die Konstanten nach der Form Vs/[(1+pT1)(1+pT2)] ein ');
    n03=input('Vs = '); T1=input('T1 = '); T2=input('T2 = ');
    d03=[(T1*T2) (T1+T2) 1];
  elseif (ss==5)|(ss==6)	
    disp(' Geben Sie die Konstanten nach der Form Vs/[(1+pT1)(1+pT2)(1+pT3] ein ');
    n03=input('Vs = '); T1=input('T1 = '); T2=input('T2 = '); T3=input('T3 = ');
    d03=[T1*T2*T3 T3*(T1+T2)+(T1*T2) T3+T2+T1 1];
  elseif (ss==7)|(ss==8)
    disp(' Geben Sie die Konstanten nach der Form Vs/(pT1) ein ');
    n03=input('Vs = '); T1=input('T1= '); d03=[T1 0];
  elseif (ss==9)|(ss==10)
    disp(' Geben Sie die Konstanten nach der Form Vs/[p(p+T1)] ein ');
    n03=input('Vs = '); T1=input('T1= '); d03=[1 T1 0];
  else
    disp(' Geben Sie die Konstanten nach der Form Vs/[p(p+T1)(p+T2)] ein ');
    n03=input('Vs = '); T1=input('T1= '); T2=input('T2= ');
    d03=[1 (T1+T2) (T1*T2) 0];
  end;
end;

% Ergebnisberechnung
% Innere Schleife
% Berechnung der Teilstrecke Gs1b
[numo1b,deno1b]=trans(n01,d01,T,m1,code1);
% innerer Regler
disp('');
disp(' Sie optimieren den inneren Regler !');
[numr1,denr1,nufil1,defil1]=habod(numo1b,deno1b);			
% Berechnung der Teilstrecke Gs1a
if co1==3    
   [numo1a,deno1a]=trans(n01,d01,T,m1,0);
elseif code1==2
   [numo1a,deno1a]=trans(n01,d01,T,m1,1);
else
   numo1a=numo1b; deno1a=deno1b;
end; end

% Zusammenfassung innere Schleife
nn=conv(numr1,numo1a); dn=conv(denr1,deno1a);
nd=conv(numr1,numo1b); dd=conv(denr1,deno1b); nd=dd+nd;
numg1=conv(nn,dd); deng1=conv(dn,nd);
[numg1f]=conv(numg1,nufil1); [deng1f]=conv(deng1,defil1);
[numg1f,deng1f]=minreal(numg1f,deng1f);

% Mittlere Schleife
% Berechnung der Teilstrecke Gs2b
nu=conv(n01,n02); de=conv(d01,d02);
[numo2b,deno2b]=trans(nu,de,T,m2,code2);
% Berechnung der offenen mittleren Schleife
nums2=conv(numg1f,numo2b); nums2=conv(nums2,deno1a);
dens2=conv(deng1f,deno2b); dens2=conv(dens2,numo1a); 
[nums2,dens2]=minreal(nums2,dens2);
e=length(nums2); f=length(dens2);
if e<f g=f-e; h=zeros(1,g); h(1,g+1)=1; nums2=conv(nums2,h);
elseif f<e g=e-f; h=zeros(1,g); h(1,g+1)=1; dens2=conv(dens2,h);
else end; end
% Mittlerer Regler
disp(' ');
disp(' Sie optimieren den mittleren Regler !');
[numr2,denr2,nufil2,defil2]=habod(nums2,dens2);	
% Berechnung der Teilstrecke Gs2a
if co2==3 
   [numo2a,deno2a]=trans(nu,de,T,m2,0);
elseif code2==2
   [numo2a,deno2a]=trans(nu,de,T,m2,1);
else
   numo2a=numo2b; deno2a=deno2b;
end; end
% Zusammenfassung mittlere Schleife
nn1=conv(numr2,numg1f); nn2=conv(nn1,numo2a);
dn1=conv(denr2,deng1f); dn2=conv(dn1,deno2a);
nd1=conv(nn1,numo2b); dd1=conv(dn1,deno2b);
nd2=conv(numo1a,dd1); nd1=conv(nd1,deno1a);
dd2=conv(dd1,deno1a);
e=length(nd1); f=length(nd2);
if e<f g=f-e; h=zeros(1,g); h(1,g+1)=1; nd1=conv(nd1,h);
elseif f<e g=e-f; h=zeros(1,g); h(1,g+1)=1; nd2=conv(nd2,h);
else end; end
nd3=nd1+nd2;
numg2=conv(nn2,dd2); deng2=conv(dn2,nd3);
[ao2,bo2,co2,do2]=tf2ss(numg2,deng2);
[w,x,y,z]=minreal(ao2,bo2,co2,do2);
[numg2,deng2]=ss2tf(w,x,y,z,1);		
[numg2f]=conv(numg2,nufil2); [deng2f]=conv(deng2,defil2);

if sch==3
  % Äußere Schleife
  % Berechnung der Teilstrecke Gs3b
  nu=conv(nu,n03); de=conv(de,d03);
  [numo3b,deno3b]=trans(nu,de,T,m3,code3);
  % Berechnung der offenen äußeren Schleife
  nums3=conv(numg2f,numo3b); nums3=conv(nums3,deno2a);
  dens3=conv(deng2f,deno3b); dens3=conv(dens3,numo2a); 
  [nums3,dens3]=minreal(nums3,dens3);
  e=length(nums3); f=length(dens3);
  if e<f g=f-e; h=zeros(1,g); h(1,g+1)=1; nums3=conv(nums3,h);
  elseif f<e g=e-f; h=zeros(1,g); h(1,g+1)=1; dens3=conv(dens3,h);
  else end
  % Äußerer Regler
  disp(' ');
  disp(' Sie optimieren den aeusseren Regler !');
  [numr3,denr3,nufil3,defil3]=habod(nums3,dens3);	
  % Berechnung der Teilstrecke Gs3a
  if co3==3 
     [numo3a,deno3a]=trans(nu,de,T,m3,0);
  elseif code3==2
     [numo3a,deno3a]=trans(nu,de,T,m3,1);
  else
     numo3a=numo3b; deno3a=deno3b;
  end
  % Zusammenfassung äußere Schleife
  nn1=conv(numr3,numg2f); nn2=conv(nn1,numo3a);
  dn1=conv(denr3,deng2f); dn2=conv(dn1,deno3a);
  nd1=conv(nn1,numo3b); dd1=conv(dn1,deno3b);
  nd2=conv(numo2a,dd1); nd1=conv(nd1,deno2a);
  dd2=conv(dd1,deno2a);
  e=length(nd1); f=length(nd2);
  if e<f g=f-e; h=zeros(1,g); h(1,g+1)=1; nd1=conv(nd1,h);
  elseif f<e g=e-f; h=zeros(1,g); h(1,g+1)=1; nd2=conv(nd2,h);
  else end
  nd3=nd1+nd2;
  numg3=conv(nn2,dd2); deng3=conv(dn2,nd3);
  [ag3,bg3,cg3,dg3]=tf2ss(numg3,deng3);
  [w,x,y,z]=minreal(ag3,bg3,cg3,dg3);
  [numg3,deng3]=ss2tf(w,x,y,z,1);		
  [numg3f]=conv(numg3,nufil3); [deng3f]=conv(deng3,defil3);
end

% graph. Ergebnisausgabe
if sch~=3  numg3f=numg2f; deng3f=deng2f; end
[a3,b3,c3,d3]=tf2ss(numg3f,deng3f) ;
[ap,bp]=d2c(a3,b3,T);
[nu,de]=ss2tf(ap,bp,c3,d3,1);
input('Wuenschen Sie eine Auswertung? nein(bel.)/ja(1): ');
disp(' ');
if ans==1 
  input('Waehlen Sie die Simulationszeit[ms] (30<sz<200) : ');
  if ans>200 disp(' Simulationszeit zu gross'),
     input('Waehlen Sie neue Zeit [ms] :');
     sz=ans;
  elseif  ans<30  disp(' Simulationszeit zu klein'),
     input('Waehlen Sie neue Zeit [ms] :');
     sz=ans;
  else sz=ans; end
  t=0:0.2:sz; n=(5*sz); u=step(nu,de,t); tt=t.'; su=0;
  for i=1:n
    if (abs(u(i+1)-u(i))>1e-4)&(u(i)>0),k=i; end
				% 1e-4 Genauigkeitsschranke
  end
  m=abs(1-u(k+1))*100; disp(' ');
  disp(' stationaere Regelabweichung in % nach   Abtastschritten :');
  disp( [m (ceil(k*.2/T))]);
  if  m>5 disp(' Reglerintegralanteil notwendig');
  elseif k==n,
    disp(' stationaerer Endzustand wurde nicht erreicht')
    disp(' -Genauigkeitsschranke vergroessern');
    disp(' -oder Berechnungszeit vergroessern');
  else [q,w]=size(u); x(1:q,1)=u; x(1:q,2)=tt;
    [hue,t1]=max(x); t11=fix(t1(1,1)/(5*T)); hue=(hue(1,1)-1)*100; 
    if hue<0.01 hue=0; t1=0; su=1;
      disp(' Die Ueberschwingweite betraegt 0 %'); disp(' ');
      disp(' Die Unterschwingweite betraegt 0 %');
    else disp(' Ueberschwingweite in % nach   Abtastschritten:');
      disp( [hue t11]); 
    end
    if su~=1
      v(1:(q-t1(1,1)),1)=x((t1(1,1)+1):q,1);
      v(1:(q-t1(1,1)),2)=x((t1(1,1)+1):q,2);
      [hun,t2]=min(v); t2=fix(t2(1,1)/(5*T)); hun=(hun(1,1)-1)*100;
      if ((abs(hun)<0.1)&(t2==1))
        hun=0; t2=0; disp(' Die Unterschwingweite betraegt 0 %');
      else disp(' Unterschwingweite in % nach   Abtastschritten :');
        disp( [hun t2+t11]); 
      end	
    end;
  end;
  input('Wuenschen Sie eine graph. Sprungantwort? nein(bel.)/ja(1): ');
  if ans==1 
    f=0:T:sz; l=(sz/T+1); v=dstep(numg3f,deng3f,l);
    plot(f,v,'ro'); grid; hold; xlabel('Zeit in ms');
    title('diskrete Sprungantwort'); k=1;
    for f=T:T:(sz-T) k=k+1; x=(v(k)-v(k-1))/50; r=T/50;
      for w=v(k-1):x:v(k) plot(f,w,'r.'); end;	
      for y=f:r:(f+T)  plot(y,v(k),'r.'); end;
    end
  pause;
  else end
else 
end;

% Ergebnisausgabe
reg1=[numr1;denr1]; fil1=[nufil1;defil1];
[g1,h1]=size(reg1); [g4,h4]=size(fil1);  disp(' ');
disp (' Der innere Reglerzaehler : '); disp( reg1(1,1:h1));
disp (' Der innere Reglernenner  : '); disp( reg1(2,1:h1));
if ((h4~=1)&(fil1(1,1)~=1))
  disp (' Der innere Vorfilterzaehler : '); disp( fil1(1,1:h4));
  disp (' Der innere Vorfilternenner  : '); disp( fil1(2,1:h4));
else end
reg2=[numr2;denr2]; fil2=[nufil2;defil2];
[g2,h2]=size(reg2); [g5,h5]=size(fil2);
disp (' Der mittlere Reglerzaehler : '); disp( reg2(1,1:h2));
disp (' Der mittlere Reglernenner  : '); disp( reg2(2,1:h2));
if ((h5~=1)&(fil2(1,1)~=1))
  disp (' Der mittlere Vorfilterzaehler : '); disp( fil2(1,1:h5));
  disp (' Der mittlere Vorfilternenner  : '); disp( fil2(2,1:h5));
else end
if sch==3
  numr3=1e10*numr3; numr3=round(numr3); numr3=numr3/1e10;
  nufil3=1e10*nufil3; nufil3=round(nufil3); nufil3=nufil3/1e10;
  defil3=1e10*defil3; defil3=round(defil3); defil3=defil3/1e10;
  reg3=[numr3;denr3]; fil3=[nufil3;defil3];
  [g3,h3]=size(reg3); [g6,h6]=size(fil3);
  disp (' Der aeussere Reglerzaehler : '); disp( reg3(1,1:h3));
  disp (' Der aeussere Reglernenner  : '); disp( reg3(2,1:h3));
  if ((h6~=1)&(fil3(1,1)~=1))
    disp (' Der aeussere Vorfilterzaehler : '); disp( fil3(1,1:h6));
    disp (' Der aeussere Vorfilternenner  : '); disp( fil3(2,1:h6));
  else end
end
