function [out1,out2,out3,out4] = eez(n,ns2,k2,Vz2)
%EEZ      Optimierung digitaler Regler auf Endliche Einstellzeit
%         ausgehend von einer Regelstrecke im Z-Bereich (trans.m)
%   ==> n/ns2: Streckenzähler bzw. -nenner im Z-Bereich
%       Es sind zwei Eingabeversionen zulässig:
%	a, nach MATLAB: [out1,out2(,out3,out4)] = eez(n,ns2)
%                               Führende Nullen in "n" und nachgestellte
%                               Nullen in "ns2" werden ausgewertet ! ns2(1)=1,
%       b, nach ROSI:   [out1,out2(,out3,out4)] = eez(n,ns2,k,Vz)
%                               Keine führenden Nullen in "n", keine nach-
%                               gestellten Nullen in "ns2" ! n(1)=ns2(1)=1,
%                               Expliziete Übergabe von Vz2 und k2
%                 Vz2           Streckenverstärkung im Z-Bereich
%                 k2            ganzzahlige Totzeit der übergebenen Strecke
%   ==> Menügesteuert kann zwischen zwei Ausgabeversionen gewählt werden:
%       a, nach MATLAB: out1 - Reglerzählerpolynom (fallende Z-Potenzen)
%                       out2 - Reglernennerpolynom (   "          "    )
%                       out3 - Filterzählerpolynom (   "          "    )
%                       out4 - Filternennerpolynom (   "          "    )
%       b, nach ROSI:   out1 - Reglerzähler bei ausgeklammertem Vr
%                       out2 - Vektor = [Vr Vf I b1 ... bn]
%                       out3 - Filterzähler mit ausgeklammertem Vf
%                       out4 - Filternenner
%                 Vr            Reglerverstärkung
%                 I             Ordnung der I-Anteile
%                 bx            Reglernennerkoeffizienten bei b0=1 und
%                               ausgeklammerten I-Anteilen
%                 Vf            Filterverstärkung
%   ==> (,out3,out4)  nur sinnvoll bei Störoptimierung

% SCHLÜSSELWORTE: Optimierung auf Endliche Einstellzeit, Digitaler Regler,
%                 Regelschleife, Störoptimierung, Führungsoptimierung,
%                 Führungsgrößenfilter, Menüführung, Kompensation wahlweise,
%                 Suboptimale Lösung, Stellgrößengleichsetzen, Ordnung des
%                 Eingangssignales

%  Copyright (c) 1997 Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%                     Tel./Fax (0351)-463-2917/7280
%                     e-mail: geitner@eeiwzb.et.tu-dresden.de
%  Literatur zur Einstellung auf Endliche Einstellzeit: Geitner, G.-H.:
%                     Entwurf digitaler Regler für elektrische Antriebe.
%                     VDE-Verlag GmbH Berlin und Offenbach 1996
%                     ISBN 3-8007-1847-2

if ((nargin<2)|(nargin==3)|(nargin>4)) error('Eingangsargumentfehler'); end
[znum,snum]=size(n); [zden,sden]=size(ns2); eps1=1e-14; % da eps zu klein ist!
if ((znum~=1)|(zden~=1))
   error('Streckenkoeffizienten als Vektoren übergeben !'); end
if ((snum>9)|(sden>9))
   disp('Ist die Ordnung der Eingangsvektoren tatsächlich so groß ?'); end
if abs(ns2(sden))<eps1 i=1; j=[1 0];
   while abs(ns2(sden-i))<eps1 i=i+1; j(i+1)=0; end; ns2=deconv(ns2,j);
   sden=sden-i; end
if abs(n(snum))<eps1 i=1; j=[1 0];
   while abs(n(snum-i))<eps1 i=i+1; j(i+1)=0; end; n=deconv(n,j);
   snum=snum-i; end
if ns2(1)~=1 error('Nenner(1)=1 erwartet!'); end; j=0;
if nargin==2
   while abs(n(j+1))<eps1; j=j+1; end;   Vz2=n(j+1);
   for i=(1+j):snum   n(i)=n(i)/Vz2;  end; k2=j;
else
   if n(1)~=1 error('Z„hler(1)=1 erwartet!'); end; end
if k2==0 error('Anwendung auf sprungfähige Strecken ist unzulässig !'); end
if (k2>=3) disp('Warnung: Ist die Totzeit tatsächlich so groß ?'); end
zs2=1; for i=(2+j):snum zs2(i-j)=n(i); end; Ns2=snum-1-j; M=sden-1;
K=[1 -1 0 0 0 0;1 -2 1 0 0 0;1 -3 3 -1 0 0;1 -4 6 -4 1 0;1 -5 10 -10 5 -1];
fos=input('Optimierung: Führung(bel.)/Störung(1) =');
  if isempty(fos)|fos~=1 fos=0;
  else disp('Voraussetzung: Störgrößeneingriff am Streckeneingang'); end;
disp('Berücksichtigung der Eingangssignalform durch Vorgabe von :');
ans=input('I-Anteilen(ET oder 1), Signalformen(2), Polynomen(3): ');
if isempty(ans) ans=1; end
if ans==1 Zqw=0; zqw=1;
   I=input('Ordnung des Reglerintegralanteils (ET=1) = ');
   if isempty(I) I=1; end
   if ((I<0)|(I>5)|(rem(I,1)~=0)) error('Keine sinnvolle Eingabe'); end
   Nw=0; nw=1; Nqw=I; nqw=1; for i=2:(I+1) nqw(i)=K(I,i); end
elseif ans==2
   I=input('Sprung(1), Rampe(2), Parabel(3): ');
   if     I==1 Zqw=0; zqw=1;
   elseif I==2 Zqw=1; zqw=[0 1];
   elseif I==3 Zqw=2; zqw=[0 0.5 0.5];
   else   error('Falsche Eingabe'); end
   Nw=0; nw=1; Nqw=I; nqw=1; for i=2:(I+1) nqw(i)=K(I,i); end
elseif ans==3
   Zqw=input('Zählerordnung= '); zqw=1; T='Eingabe von z0= ';
   for i=1:Zqw T(14)=48+i; T=setstr(T); zqw(1+i)=input(T); end
   Nqw=input('Nennerordnung= '); nqw=1; T='Eingabe von n0= ';
   for i=1:Nqw T(14)=48+i; T=setstr(T); nqw(1+i)=input(T); end
   I=0; j=0; hiv=nqw; while j==0  nw=hiv; [hiv,r]=deconv(hiv,[1 -1])
                            if j==sum(r) I=I+1; end; end; Nw=Nqw-I;
else   error('Falsche Eingabe'); end
pz=roots(ns2); Is=0;
for i=1:M
    if pz(i)==1 Is=Is+1; end; end
if Is>0 disp('Die Anzahl der Strecken-I-Anteile beträgt Is='); disp(Is);
   if I>0
      ans=input('Verwendung zum Reglerentwurf? nein(0)/ja(Anzahl)');
      if ans~=0
         if ((rem(ans,1)~=0)|(ans>Is)|(ans<0))
            error('Keine sinnvolle Eingabe'); end; hiv=I;I=I-ans;
         if I<0  I=0;
            disp('Die Anzahl der berücksichtigten I-Anteile beträgt:');
            disp(hiv); end
         if fos==1
            disp('Warnung: Das Ergebnis ist eine Näherungslösung !!!'); end
end; end; end
if fos==0
ans=input('"DGG" zwischen Regler und kontinuierlicher Strecke einfügen? nein(bel.)/ja(1): ');
else
ans=input('"DGG" zwischen Regler und Störgrößeneingriff einfügen ? nein(bel.)/ja(1): ');
end
if isempty(ans)|ans~=1 ans=0; end
if ans==1
disp('Definition des Differenzengleichungsgliedes: Z1 = Vz1*(1+g1*z~1+...)*z~k1');
   Vz1=input('Eingabe der Verstärkung Vz1= ');
   k1 =input('Eingabe der Totzeit k1= ');
   Ns1=input('Eingabe der Polynomordnung= ');
   if (isempty(Vz1)|isempty(k1)|isempty(Ns1)|(rem(Ns1,1)~=0))
                                             error('Falsche Eingabe'); end
   T='Eingabe von g0= '; zs1=1;
   for i=1:Ns1 T(14)=48+i; T=setstr(T); zs1(1+i)=input(T); end
else Vz1=1; k1=0; Ns1=0; zs1=1; end; k=k1+k2; Vz=Vz1*Vz2;
if fos==1    % Beginn Störoptimierung
Ns2str=0; zs2str=1;
ans=input('Nullstellenkompensation: ja(bel.)/nein(1) =');
  if isempty(ans)|ans~=1 ans=0; end
  if ans==1 Ns2ste=Ns2; zs2ste=zs2;
  else
    if Ns2==0 error('Diese Strecke besitzt keine Nullstellen ~=0'); end
    nz=roots(zs2); l=Ns2;
    disp('Die Strecke besitzt folgende Nullstellen:'); disp(nz)
    Ns2ste=0; zs2ste=1;
    while l~=0
       if imag(nz(l))==0
          disp('Kompensation der Nullstelle: '); disp(nz(l))
          anzn=1;
       else
          disp('Kompensation der Nullstelle: '); disp(nz(l))
                                                 disp(nz(l-1))
          anzn=2;
       end
       ans=input('ja(bel.)/nein(1)= ');
       if isempty(ans)|ans~=1 ans=0; end
       for i=1:anzn
           if ans~=1
              Ns2str=Ns2str+1; zs2str(Ns2str+1)=0; arbvek=1; for j=1:Ns2str
                               arbvek(j+1)=zs2str(j+1)-zs2str(j)*nz(l+1-i);
                               end; zs2str=arbvek;
           else
              Ns2ste=Ns2ste+1; zs2ste(Ns2ste+1)=0; arbvek=1; for j=1:Ns2ste
                               arbvek(j+1)=zs2ste(j+1)-zs2ste(j)*nz(l+1-i);
                               end; zs2ste=arbvek;  end
       end; l=l-anzn; end   %"while" ist sinnvoll da komplexe Nullstellen
                            %im Paar kompensiert werden sollen
  if Ns2ste~=0
  disp('Hinweis: Minimale Reglerordnung nur durch Vollkompensation erreichbar');
  end; end
mst=k1+2*k2+Ns1+Ns2+Ns2ste+Zqw; ast=1;
disp('Die minimale Stör-Anregelschrittzahl beträgt: '); disp(mst);
Ast=input('Anzahl zusätzlicher Stör-Anregelschritte (subopt. Einstellung, ET=>0) =');
 if isempty(Ast) Ast=0; end
 if ((Ast<0)|(Ast>4)|(rem(Ast,1)~=0)) error('Keine sinnvolle Eingabe'); end
 if Ast~=0 ans=input('Vorgabe der As-Polynomkoeffizienten: maschinell(bel.)/Eingabe(1):');
    if isempty(ans)|ans~=1 ans=0; end
    if ans==1
       T='Eingabe von a0 =';
       for i=1:Ast T(14)=48+i; T=setstr(T); ast(i+1)=input(T); end
    else
     yg=ones(1,Ast+1); Er=Ns2ste+Zqw-Nqw;
     if Er>0 yg(Ast+1+Er)=0; end
     a1=deconv(conv(yg,nqw),conv(zs2ste,zqw));
    disp('Hinweis: Die Berechnng erfolgt durch Stellgrößengleichsetzen')
       for i=1:Ast ast(1+i)=a1(1+i); end; end; end % wegen Polynomdivision nötig
AV=Nqw-1+M+Ast; BV=k+Ns1+Ns2ste-1+Ast;
cp=conv(nqw,ns2);     CP=Nqw+M;      % mit cp0=1, CP echte Ordnung
dp=conv(zs1,zs2ste);  DP=Ns1+Ns2ste; %  "  dp0=1, DP   "      "
if BV>0
   hiv=-1; cpe=-cp;
     if BV>1  hiv(BV)=0; cpe(BV+CP)=0; end
   ma=toeplitz(cpe,hiv); end
hiv=-conv(dp,ast); HIV=DP+Ast;        % mit hiv0=1
for i=0:HIV
    ma(k+i,BV+1)=hiv(i+1); end
for j=0:(AV-Ast-1)
    for i=0:DP
        ma(k+i+j+1+Ast,BV+2+j)=-dp(i+1); end; end
hiv=0; hiv(BV+CP)=0;
for i=1:CP
    hiv(i)=cp(i+1); end
erg=inv(ma)*hiv';
bv=1; if BV>0
         for i=1:BV
             bv(1+i)=erg(i); end; end
av=erg(1+BV);                        % entspricht V0
for i=1:Ast
    av(1+i)=ast(1+i)*av(1); end
for i=1:(AV-Ast)
    av(1+Ast+i)=erg(1+BV+i); end
out1=av/Vz;                          % Reglerzähler mit einmultipliziertem Vr
out2=conv(conv(bv,zs2str),nw);       % I-Anteile werden später zugerechnet
ver=BV+Ns2str+Nw;
end
% Führungsoptimierung
Ms2str=0; ns2str=1;
if fos==0 ans=input('Führungsoptimierung mit Polkompensation: ja(bel.)/nein(1)');
else      ans=input('Filterentwurf mit Polkompensation: ja(bel.)/nein(1) ='); end
  if isempty(ans)|ans~=1 ans=0; end
  if ans==1 Ms2ste=M; ns2ste=ns2;
  else
    if M==0 error('Diese Strecke besitzt keine Pole'); end
    pz=roots(ns2); l=M;
    disp('Die Strecke besitzt folgende Pole:'); disp(pz)
    Ms2ste=0; ns2ste=1;
    while l~=0
       if imag(pz(l))==0
          disp('Kompensation des Poles: '); disp(pz(l))
          anzn=1;
       else
          disp('Kompensation des Poles: '); disp(pz(l))
                                            disp(pz(l-1))
          anzn=2;
       end
       ans=input('ja(bel.)/nein(1)= ');
       if isempty(ans)|ans~=1 ans=0; end
       for i=1:anzn
           if ans~=1
              Ms2str=Ms2str+1; ns2str(Ms2str+1)=0; arbvek=1; for j=1:Ms2str
                               arbvek(j+1)=ns2str(j+1)-ns2str(j)*pz(l+1-i);
                               end; ns2str=arbvek;
           else
              Ms2ste=Ms2ste+1; ns2ste(Ms2ste+1)=0; arbvek=1; for j=1:Ms2ste
                               arbvek(j+1)=ns2ste(j+1)-ns2ste(j)*pz(l+1-i);
                               end; ns2ste=arbvek;  end
       end; l=l-anzn; end   %"while" ist sinnvoll da komplexe Nullstellen
                            %im Paar kompensiert werden sollen
  if Ms2ste~=0
  if fos==0
  disp('Hinweis: Minimale Reglerordnung nur durch Vollkompensation erreichbar');
  else
  disp('Hinweis: Minimale Filterordnung nur durch Vollkompensation erreichbar');
  end; end; end
mfu=k+Ns1+Ns2+Ms2ste+Zqw; afu=1;
disp('Die minimale Führungs-Anregelschrittzahl beträgt: '); disp(mfu);
Afu=input('Anzahl zusätzlicher Führungs-Anregelschritte (subopt. Einstellung, ET=>0) =');
 if isempty(Afu) Afu=0; end
 if ((Afu<0)|(Afu>4)|(rem(Afu,1)~=0)) error('Keine sinnvolle Eingabe'); end
 if Afu~=0 ans=input('Vorgabe der Af-Polynomkoeffizienten: maschinell(bel.)/Eingabe(1):');
    if isempty(ans)|ans~=1 ans=0; end
    if ans==1
       T='Eingabe von a0 =';
       for i=1:Afu T(14)=48+i; T=setstr(T); afu(i+1)=input(T); end
    else
     yg=ones(1,Afu+1); Er=M+Zqw-Nqw;
     if Er>0 yg(Afu+1+Er)=0; end
     a1=deconv(conv(yg,nqw),conv(ns2,zqw));
    disp('Hinweis: Die Berechnng erfolgt durch Stellgrößengleichsetzen')
       for i=1:Afu afu(1+i)=a1(1+i); end; end; end % wegen Polynomdivision nötig
AVf=Nqw-1+Ms2ste+Afu; BVf=k+Ns1+Ns2-1+Afu;
cp=conv(nqw,ns2ste);  CP=Nqw+Ms2ste;      % mit cp0=1, CP echte Ordnung
dp=conv(zs1,zs2);     DP=Ns1+Ns2;         %  "  dp0=1, DP   "      "
if BVf>0
   hiv=-1; cpe=-cp;
      if BVf>1  hiv(BVf)=0; cpe(BVf+CP)=0; end
   ma=toeplitz(cpe,hiv); end
hiv=-conv(dp,afu); HIV=DP+Afu;        % mit hiv0=1
for i=0:HIV
    ma(k+i,BVf+1)=hiv(i+1); end
for j=0:(AVf-Afu-1)
    for i=0:DP
        ma(k+i+j+1+Afu,BVf+2+j)=-dp(i+1); end; end
hiv=0; hiv(BVf+CP)=0;
for i=1:CP
    hiv(i)=cp(i+1); end
erg=inv(ma)*hiv';
bvf=1; if BVf>0
         for i=1:BVf
             bvf(1+i)=erg(i); end; end
avf=erg(1+BVf);                        % entspricht V0
for i=1:Afu
    avf(1+i)=afu(1+i)*avf(1); end
for i=1:(AVf-Afu)
    avf(1+Afu+i)=erg(1+BVf+i); end
out1f=conv(ns2str,avf/Vz);           % Reglerzähler mit einmultipliziertem Vr
out2f=conv(bvf,nw);                  % I-Anteile werden später zugerechnet
verf=BVf+Nw;
if fos==1 % Filterberechnung
   out3=conv(zs2str,avf)/av(1);
   out4=out1/out1(1);
else out3=1; out4=1; out1=out1f; out2=out2f; ver=verf;
end
ans=input('Ergebnisausgabe nach MATLAB(bel.)/ROSI(1): ');
if isempty(ans)|ans~=1 ans=0; end
if ans==1 out2=conv(out2,nw); ver=ver+Nw;
          for i=1:ver out2(ver+4-i)=out2(ver+2-i); end
          out2(1)=out1(1); out2(2)=out4(1); out2(3)=I;
          out1=out1/out1(1);
          out4=out4/out4(1);
else      out2=conv(out2,nqw);
end
