function [out1,out2,out3,out4] = bod_2(inp1,inp2,inp3)
%***** same function as bod.m, BUT application of OLD fsolve calling syntax *****

%BOD      Berechnung eines digitalen Reglers nach dem Digitalen Betragsoptimum
%         (BOD) - ausgehend von einer Regelstrecke im Z-Bereich (trans.m).
%         LTI Syntax möglich; Menüführung abschaltbar.
%   ==> Es sind zwei Eingabeversionen zulässig:
%       a, alte MATLAB Syntax: Zähler- und Nennerpolynom
%          [out1,out2(,out3,out4)] = bod_2(inp1,inp2(,inp3))
%                       inp1 - Zähler der Z-transformierten Strecke
%                       inp2 - Nenner der Z-transformierten Strecke
%          Führende Nullen in "inp1" und nachgestellte Nullen in "inp2" werden 
%          ausgewertet! Nenner(1)=1 wird erwartet.
%       b, neue MATLAB Syntax: LTI-Modell
%          [out1(,out2)] = bod_2(inp1(,inp2))
%                       inp1 - LTI-Modell der Z-transformierten SISO Strecke
%   ==> Entsprechend der Syntax erfolgt die Ergebnisausgabe in zwei Versionen:
%       a, alte Syntax: out1 - Reglerzählerpolynom (fallende Z-Potenzen)
%                       out2 - Reglernennerpolynom (    "         "    )
%                       out3 - Filterzähler
%                       out4 - Filternennerpolynom (    "         "    )
%       b, neue Syntax: out1 - LTI-Modell Regler
%                       out2 - LTI-Modell Filter
%   ==> (,out3,out4) bzw. (,out2) ist nur bei Störoptimierung sinnvoll.
%   ==> inpx = inp3(alte Syntax) bzw. inp2(LTI) = Abschaltung der Menüführung:
%       inpx(1) - bel./1 = Optimierung: Führungsverhalten/Störverhalten
%       inpx(2) - bel./1 = Polkompensation: nein/ja
%       inpx(3) - bel./1 = Reglerintegralanteil: ja/nein
%       inpx(4) = Anzahl zu kompensierender Pole (komplex: Polpaar; Polreihenfolge
%                 Betrag: abnehmend bei Fürungs- / zunehmend bei Störoptimierung)
%       inpx(5) = Ordnung des zu >>berechnenden<< Reglerzähleranteils
%                 (Gesamtzählerordnung: Kompensation+Berechnung)
%       inpx(6) - bel./1/2 = Unterdrückung von Ausgabetexten: nein/ja(K)/ja(W+K)
%                 (W - Warnung; K - Kommentar)
%       inpx(7) - bel./1 = Sonderfall gedämpfter D-Anteil einschließen: nein/ja
%       inpx(8) - 0/1/2 - Auswahl aus drei Lösungen der kub. Gl. für b1 (inpx(7)=1)
%       inpx(9) - bel./1 =Lösung über quadratische Gleichung - Auswahl: V02/V01
%                 (Erfahrungswert: V02 als Verstärkung des offenen Kreises)
%       inpx(10)- bel./1/2 =Lösung über Suchalgorithmus - Vorgabe von Anfangswerten:
%                 maschinell/von Hand PID-Reglerparameter/von Hand Reglernullstelle
%       inpx(11)...inpx(13)=ohne Bedeutung   für inpx(10)=bel.
%                          =[VR d1 d2]        "  inpx(10)=1
%                          =[Nullstelle % %]  "  inpx(10)=2
%   bel. - beliebiger anderer Wert, z.B. 0, entspricht Vorzugseinstellung;
%   >>>Die Abschaltung von Ausgabetexten wird erst nach Einarbeitung empfohlen!<<<
%   Default und Ergänzung fehlender Elemente von inpx: Null

% SCHLÜSSELWORTE: Digitales Betragsoptimum, Digitaler Regler von PID-Typ,
%		            Regelschleife, Störoptimierung, Führungsoptimierung,
%		            Führungsgrößenfilter, Symmetrisches Optimum, Menüführung wahlweise,
%		            Polkompensation wahlweise, LTI

%  ***** Voraussetzungen: Vorhandensein von bodg, vobo ****
%  Copyright (c) 1997-2007. All Rights Reserved.
%                     Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%                     http://eeiwzg.et.tu-dresden.de/ae2_files/ae_1.htm
%  Literatur zum Digitalen Betragsoptimum: Geitner, G.-H.:
%                     Entwurf digitaler Regler für elektrische Antriebe.
%                     VDE-Verlag GmbH Berlin und Offenbach 1996
%                     ISBN 3-8007-1847-2

FT01='Eingangsargumentefehler';
FT02='Alte Syntax: Streckenkoeffizienten als Vektoren übergeben';
FT03='Nenner(1)=1 erwartet!';
FT04='Eingangsargument Steuervektor(inpx) - falsche Dimension';
FT05='Eingangsargument Steuervektor(inpx) - Element 4 oder/und 5 fehlerhaft!';
FT06='Polkompensation unmöglich: Die eingegebene Strecke besitzt keine Pole';
FT07='Mit dieser Einstellung ist keine Berechnung möglich!';
FT08='Nicht sinnvoll - zwei Freiheitsgrade zur Störoptimierung notwendig';
FT09='Falsche Eingabe zur Auswahl einer Lösung für b1';
FT10='Auswahl von Lösungen für b1: Komplexe Reglerparameter sind unzulässig';
FT11='Störoptimierung nicht möglich - zu niedrige Streckenstruktur';
FT12='Störopt. PID-Regler: nicht genügend Pole für maschinelle Anf.-werte';
FT13='Störoptimierung: Reglergrad > PID ist nicht vorbereitet!';
WT01='Warnung: Ist die Ordnung der Eingangsvektoren tatsächlich > 9 ?';
WT02='Warnung: Algorithmus wird auf sprungfähige Strecke angewendet';
WT03='Warnung: Ist die Totzeit tatsächlich > 3 ?';
WT04='Warnung: Mit dieser Einstellung ist keine optimale Lösung möglich!';
WT05='Warnung: Polanzahl geringer als Anzahl zu kompensierender Pole!';
WT06='Warnung: Automatische Erhöhung der Reglerordnung wegen komplexer Pole!';
WT07='Warnung: Vorhandensein eines I-Anteiles vor Störmischstelle Überprüfen!';
WT08='Warnung: Mit gedämpftem D-Anteil - Nur in Sonderfällen brauchbar!!';
WT09='Warnung: Ist diese hohe Reglerordnung > 4 sinnvoll?';
WT10='Warnung: nur  -1 < b1 < 0  entspricht einer Dämpfungszeitkonstante';
KT01='Die Strecke besitzt folgende Pole:';
KT02='Kompensation des Poles: ';
KT03='Kompensation der Pole: ';
KT04='Die kubische Gleichung für b1 liefert folgende Lösungen:';
KT05='Störoptimierung: Lösungen für Verstärkung des offenen Kreises: (V01 V02)';
KT06='Störoptimaler PID-Regler: Lösungsalgorithmus benötigt Anfangswerte!';
KT07='Anfangswert für die Verstärkung des offenen Kreises: (V01 V02)';
KT08='Anfangswerte des Suchalgorithmus (VR d1 d2) :';

if (nargin<1)|(nargin>3)|(nargin==2&isa(inp2,'lti'))|...
   (nargin==3&isa(inp3,'lti')) error(FT01); end
if isa(inp1,'lti')          %LTI-Modell übergeben
   if (nargin==3)|~issiso(inp1)|~isdt(inp1) error(FT01); end
   if ~isa(inp1,'tf') inp1=tf(inp1); end
   [n,d]=tfdata(inp1,'v'); LTI=1;
   if nargin==2 MUEF=0; STEUV=inp2; else MUEF=1; STEUV(6)=0; end   
else                        %alte Syntax übergeben
   if nargin==1 error(FT01); end
   n=inp1; d=inp2; LTI=0;
   if nargin==3 MUEF=0; STEUV=inp3; else MUEF=1; STEUV(6)=0; end
end
[znum,snum]=size(STEUV); if znum>1|snum>13 error(FT04); end
if snum<13 STEUV(13)=0; end %Default-Ergänzung
if STEUV(4)<0|rem(STEUV(4),1)~=0 error(FT05); end
if STEUV(5)<0|rem(STEUV(5),1)~=0 error(FT05); end
          %für restliche STEUV-Prüfung kann Eingabeprüfung verwendet werden!
[znum,snum]=size(n); [zden,sden]=size(d); eps1=1e-14; % da eps zu klein ist!
if ((znum~=1)|(zden~=1)) error(FT02);  end;
if STEUV(6)~=2 if ((snum>9)|(sden>9)) disp(WT01); end; end
if abs(d(sden))<eps1 i=1; j=[1 0];
   while abs(d(sden-i))<eps1 i=i+1; j(i+1)=0; end; d=deconv(d,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 d(1)~=1 error(FT03); end
j=0; while abs(n(j+1))<eps1 j=j+1; end; Vz=n(j+1);
     for i=(1+j):snum   n(i)=n(i)/Vz; end; k=j;
if STEUV(6)~=2 if k==0 disp(WT02); end; if k>3 disp(WT03); end; end

for i=(2+j):snum  siz(1,1+i-j)=n(i); end
siz(1,1)=Vz; siz(1,2)=snum-1-j;
if   MUEF~=0 ans=input('Optimierung: Führung(bel.)/Störung(1) ='); 
else ans=STEUV(1); end
  if isempty(ans)|ans~=1 ans=0; end 
  opa(1)=ans;
if   MUEF~=0 ans=input('Polkompensation: nein(bel.)/ja(1) =');
else ans=STEUV(2); end    
  if isempty(ans)|ans~=1 ans=0; end
  opa(2)=ans; DKo=0; RZKo=1;
  if ans==1
     if opa(1)==1&STEUV(6)~=2 disp(WT04); end %Störoptimierung
  if sden==1 error(FT06); end
    pz=roots(d);
    if opa(1)==0 pz=sort(pz); end %STEUV: Unterschied stör- o. fürungsopt.
    l=sden-1;
    if STEUV(6)~=2 if l<STEUV(4) disp(WT05); end; end
    if STEUV(6)~=2&STEUV(6)~=1 disp(KT01); disp(pz); end
    hid=1; shid=1; sRZKo=1; % RZKo=1 ist schon gesetzt
    while l~=0 %Whileschleife zur Polkompensation
       if imag(pz(l))==0
          if STEUV(6)~=2&STEUV(6)~=1 disp(KT02); disp(pz(l)); end; anzp=1;
       else
          if STEUV(6)~=2&STEUV(6)~=1 disp(KT03); disp(pz(l));
                                         disp(pz(l-1)); end; anzp=2;
       end
       if   MUEF~=0 ans=input('nein(bel.)/ja(1)= ');
            if isempty(ans)|ans~=1 ans=0; end
       else ans=STEUV(4);
            if ans~=0 ans=ans-anzp; 
               if ans<0 STEUV(4)=0; if STEUV(6)~=2 disp(WT06); end
               else  STEUV(4)=ans;
               end; ans=1; end
       end
       if MUEF==0&STEUV(6)~=2&STEUV(6)~=1 
          if ans==0 disp('Nein'); else disp('Ja'); end; end
       for i=1:anzp
           if ans==1
              sRZKo=sRZKo+1; RZKo(sRZKo)=0; arbvek=1; for j=1:(sRZKo-1)
                               arbvek(j+1)=RZKo(j+1)-RZKo(j)*pz(l+1-i);
                               end; RZKo=arbvek;
           else
              shid=shid+1; hid(shid)=0; arbvek=1; for j=1:(shid-1)
                               arbvek(j+1)=hid(j+1)-hid(j)*pz(l+1-i);
                               end; hid=arbvek;  end
       end; l=l-anzp; 
     end %"while" ist sinnvoll da komplexe Pole im Paar kompensieren!
  sden=shid; d=hid; DKo=sRZKo-1; end %RZKo liegt berechnet vor
for i=2:sden    siz(2,1+i)=d(i);  end
siz(2,1)=k; siz(2,2)=sden-1;
if   MUEF~=0 ans=input('Reglerintegralanteil: ja(bel.)/nein(1) =');
else ans=STEUV(3); end
  if isempty(ans)|ans~=1 ans=1; else ans=0; end
  opa(4)=ans;
  if ((ans==0)&(opa(1)==1)) if STEUV(6)~=2 disp(WT07); end; end
  if ((ans==0)&(siz(2,2)==0)) error(FT07); end %Keine Berechnung möglich 
if   MUEF~=0 ans=input('Ordnung des zu berechnenden Reglerzähleranteils (ET=0) =');
             if isempty(ans) ans=0; end
else ans=STEUV(5); end
  D=ans; b1=0;
  if ((ans==0)&(opa(1)==1)) error(FT08); end %zwei Freiheitsgrade für Störopt. nötig
  if((ans==0)&(opa(1)==0)&(DKo>1))
    if MUEF~=0 ans=input('"gedämpfter" D-Anteil: nein(bel.)/ja(1) =');
    else ans=STEUV(7); end
    if isempty(ans)|ans~=1 ans=0;
    else if STEUV(6)~=2 disp(WT08); end; end %Nur in Sonderfällen brauchbar
    opa(5)=ans;
  else opa(5)=0; end
opa(3)=D; RZO=DKo+D;
if RZO>4 disp(WT09); end %Hohe Reglerordnung sinnvoll?
[R]=vobo(siz,opa);
if opa(1)==0             %Führungsoptimierung
   if opa(5)==0          %ohne Dämpfung des D-Anteiles
      for i=1:(D+1)
          svek(i,1)=-R(i,1);
          for j=1:(D+1)
              RM(i,j)=R(i,j+1);
               end; end
      erg=(inv(RM)*svek)'/siz(1,1);
   else %Mit Dämpfung des D-Anteiles - Zu berechnender RZ-Anteil D=0 vorgegeben!
   K0=R(1,1)*R(2,2)-R(1,2)*R(2,1);
   K1=R(1,3)*R(2,2)-R(1,2)*R(2,3);
   K2=R(1,3)*R(2,4)-R(1,4)*R(2,3);
   K3=R(1,1)*R(2,4)-R(2,1)*R(1,4);
   pz=roots([1 (K2+K0)/K3 1+(K1/K3) K0/K3]); pz=sort(pz);
   if STEUV(6)~=2&STEUV(6)~=1 disp(KT04); disp(pz); end
   if MUEF~=0 ans=input('Auswahl einer Lösung (1,2,3): ');
   else ans=STEUV(8)+1; end
   if isempty(ans) ans=1; end          %versehentliches ET als '1' abblocken
   if (ans<1)|(ans>3) error(FT09)      %Falsche Eingabe
   elseif imag(pz(ans)~=0) error(FT10) %Komplexe Reglerparameter unzulässig
   elseif (pz(ans)<=-1)|(pz(ans)>0) disp(WT10); end; b1=pz(ans);
   erg=-(R(2,1)+b1*(R(2,3)+b1*R(2,1)))/((R(2,2)+b1*R(2,4))*siz(1,1));
   end
else   % Störoptimierung
   if D==1 %störoptimaler PI-Regler
      pzae=(R(1,4)-R(1,5))*R(2,1)-(R(2,4)-R(2,5))*R(1,1)+R(1,2)*R(2,3)...
                                                           -R(1,3)*R(2,2);
      pnen=(R(1,4)-R(1,5))*R(2,2)-(R(2,4)-R(2,5))*R(1,2);
      if pnen==0 error(FT11); end       %Streckenordnung zuniedrig
      pe=pzae/pnen; qu=(R(1,1)*R(2,3)-R(2,1)*R(1,3))/pnen;
      V0=roots([1 pe qu])';
      if STEUV(6)~=2&STEUV(6)~=1 disp(KT05); disp(V0); end %Lösungen für offenen Kreis
      if MUEF~=0 ans=input('Auswahl: V02(bel.)/V01(1)= ');
      else ans=STEUV(9); end                       %Erfahrungswert: V01
      if isempty(ans)|ans~=1 erg=V0(2); else erg=V0(1); end
      erg(2)=(erg(1)*R(2,2)+R(2,1))/(erg(1)*(R(2,4)-R(2,5))-R(2,3));
      erg=erg/siz(1,1);
   elseif D==2 %störoptimaler PID-Regler
      if STEUV(6)~=2&STEUV(6)~=1 disp(KT06); end
      if MUEF~=0 ans=input('Vorgabe: maschinell(bel.)/von Hand(1) =');
      else ans=STEUV(10); if ans==2 ans=1; end; end %Zusammenfassung für MUEF=0
      if isempty(ans)|ans~=1 ans=0; end
      if ans==1
if MUEF~=0 
   ans=input('Vorgabe einer: PID-Reglereinstellung(bel.)/Reglernullstelle(1) =');
else ans=STEUV(10)-1; end
if isempty(ans)|ans~=1 ans=0; end
        if ans==1
           if MUEF~=0 ans=input('Reglernullstelle ='); else ans=STEUV(11); end
           trz=[1 ans]; j=1; n=conv(n,trz); siz(1,2)=snum;
           for i=3:(snum+2)  siz(1,i)=n(i-1);  end
        else
           if MUEF~=0 ans=input('Anfangswert für VR= '); else ans=STEUV(11); end
           anfv(1)=ans*siz(1,1);
           if MUEF~=0 ans=input('Anfangswert für d1= '); else ans=STEUV(12); end
           anfv(2)=ans*anfv(1);
           if MUEF~=0 ans=input('Anfangswert für d2= '); else ans=STEUV(13); end
           anfv(3)=ans*anfv(1); j=0;
        end
      else    spz=size(d)-1;
         if spz==0  error(FT12); end  %Strecke besitzt keine Pole (mehr)!
         pz=roots(d); kpl=sum(abs(imag(pz)));
         if kpl==0  i1=1;
            for i=1:spz
                if pz(i)>=0
                   pzpos(i1)=pz(i); i1=i1+1;
                end; end
            if     i1>1     trz=[1 -min(pzpos)];
            else   trz=[1 -max(pz)];  end;  j=1;
            d=deconv(d,trz); siz(2,2)=sden-1-j; d(sden)=0;
            for i=3:(sden+1) siz(2,i)=d(i-1); end
         else
            anfv=[siz(1,1) siz(1,1) siz(1,1)]; j=0;
         end
      end
      opa(3)=1; [E]=vobo(siz,opa);
      if j==1
        pzae=(E(1,4)-E(1,5))*E(2,1)-(E(2,4)-E(2,5))*E(1,1)+E(1,2)*E(2,3)...
                                                              -E(1,3)*E(2,2);
        pnen=(E(1,4)-E(1,5))*E(2,2)-(E(2,4)-E(2,5))*E(1,2);
        pe=pzae/pnen; qu=(E(1,1)*E(2,3)-E(2,1)*E(1,3))/pnen;
        V0=roots([1 pe qu]);
        if STEUV(6)~=2&STEUV(6)~=1 disp(KT07); disp(V0); end %AW 'V' offener Kreis
        if MUEF~=0 ans=input('Auswahl: V02(bel.)/V01(1) =');
           else ans=STEUV(9); end
        if isempty(ans)|ans~=1 erg=V0(2); else erg=V0(1); end
        erg(2)=(erg(1)*E(2,2)+E(2,1))/(erg(1)*(E(2,4)-E(2,5))-E(2,3));
        anfv=conv(erg,trz);
        if STEUV(6)~=2&STEUV(6)~=1 disp(KT08); %AW=[VR d1 d2] des Suchalgorithmus
           disp([anfv(1)/siz(1,1) anfv(2)/anfv(1) anfv(3)/anfv(1)]); end
      end
      erg=fsolve('bodg',anfv',[0],[],R);  erg=erg'/siz(1,1);
   else error(FT13);               % "D" wurde nicht mit "1" oder "2" vorgegeben
   end; rzs=sum(erg);
end
out1=conv(erg,RZKo);                                    %Reglerzähler
if opa(1)==1 [j,i1]=size(erg);                          %Störoptimierung
             out3(i1)=0; out3(1)=rzs/out1(1);           %Filterzähler
             out4=erg/out1(1);                          %Filternenner
else out3=1; out4=1; end                                %Führungsoptimierung
if    ((-opa(4)*b1)~=0)  out2=[1 b1-opa(4) -opa(4)*b1]; %Reglernenner
elseif ((b1-opa(4))~=0)  out2=[1 b1-opa(4)];
else                     out2=[1];  end
if LTI==1 out1=tf(out1,out2,inp1,'Variable','z^-1'); %Eigenschaften wie Strecken-
          out2=tf(out3,out4,inp1,'Variable','z^-1'); %LTI-Darstellung, z.B.: T
          out3=[]; out4=[];
else      [j,i1]=size(out1); [j,i2]=size(out2);
          if i1>i2 out2(i1)=0; elseif i1<i2 out1(i2)=0; end
end
