Code covered by the BSD License

# Toolbox BOD Version 2.5

23 Mar 2001 (Updated 12 Jun 2012)

Digital Amplitude Optimum (BOD) for discontinuous control

bod.m
```function [out1,out2,out3,out4] = bod(inp1,inp2,inp3)
%BOD      Calculation of a discontinuous controller using the
%         Discontinuous Amount Optimum (BOD) - starting via a plant
%         given in z-domain, typically applying function trans of the BOD toolbox.
%         LTI syntax optionally; menu driven optionally;
%   ==> There are two input versions valid:
%       a, standard MATLAB input syntax: numerator and denominator polynomial
%          [out1,out2(,out3,out4)] = bod(inp1,inp2(,inp3))
%                       inp1 - numerator   of the z-transformed plant
%                       inp2 - denominator of the z-transformed plant
%          Leading zeros of "inp1" and zeros at the end of "inp2" are evaluated!
%          Denominator(1)=1 is expected!
%       b, LTI model MATLAB syntax:
%          [out1(,out2)] = bod(inp1(,inp2))
%                       inp1 - LTI model of the z-transformed SISO plant
%   ==> According input syntax there are two versions for the output of results:
%       a, standard syntax:
%                out1 - controller numerator   polynomial (descending powers of z)
%                out2 - controller denominator polynomial (    "         "   "  ")
%                out3 - pre-filter numerator   polynomial (    "         "   "  ")
%                out4 - pre-filter denominator polynomial (    "         "   "  ")
%       b, LTI syntax: out1 - LTI model - controller
%                      out2 - LTI model - pre-filter
%   ==> (,out3,out4) resp. (,out2) - for delayed input signals useful only.
%   ==> inpx = inp3(standard syntax) resp. inp2(LTI syntax) = menu switch off via
%                                                         following control vector:
%       inpx(1) - any/1 = Optimization regarding: undelayed / delayed input signals
%       inpx(2) - any/1 = pole compensation: no/yes
%       inpx(3) - any/1 = controller integral part: yes/no
%       inpx(4) = number of poles to be compensated (complex: pole pairs;
%                                                    real: poles - and sequence via
%                 absolute value: descending undelayed / ascending delayed -inputs)
%       inpx(5) = order of the controller numerator part to be >>calculated<<
%                 (total numerator order: compensation + calculation)
%       inpx(6) - any/1/2 = elimination of text outputs: no / yes(C) / yes(W+C)
%                 (W - warnings; C - comments)
%       inpx(7) - any/1 = implicate special case damped derivative part: no/yes
%       inpx(8) - 0/1/2 - choice from 3 cubic equa. solutions as to b1 (inpx(7)=1)
%       inpx(9) - any/1 = solution via quadratic equation - choice: V02/V01
%                 (rule of thumb: V02 as the gain of the open loop)
%       inpx(10)- any/1/2 = solution via search algorithm - initial values for PID:
%                 automatic / manual controller parameters / manual controller zero
%       inpx(11)...inpx(13) = without meaning if inpx(10)=any
%                           = [VR d1 d2]      if inpx(10)=1
%                           = [zero % %]      if inpx(10)=2
%   any == any other input, e.g. enter or 0, matches preferential input setting;
%   >>The elimination of text outputs is recommended after getting familiar only!<<
%   Defaults and completions of missing elements of inpx: zero

% KEY WORDS: Discontinuous Amplitude Optimum, discontinuous controller of
%            PID-type family, control loop, optimization regarding undelayed
%		     input signals, optimization regarding delayed input signals,
%            pre-filter design, optionally pole compensation,
%		     optionally menu driven, optionally LTI syntax

%  ***** Preconditions: existence of bodg, vobo - see BOD toolbox ****
%                     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
%  Reference regarding the application of the Discontinuous Amplitude Optimum:
%                     Geitner, G.-H.:
%                     Entwurf digitaler Regler fr elektrische Antriebe.
%                     publisher: VDE-Verlag GmbH Berlin und Offenbach 1996
%                     ISBN 3-8007-1847-2

%----------------------------German-------------------------------------------
%BOD      Berechnung eines digitalen Reglers nach dem Digitalen Betragsoptimum
%         (BOD) - ausgehend von einer Regelstrecke im Z-Bereich (trans.m).
%         LTI Syntax mglich; Menfhrung abschaltbar.
%   ==> Es sind zwei Eingabeversionen zulssig:
%       a, alte MATLAB Syntax: Zhler- und Nennerpolynom
%          [out1,out2(,out3,out4)] = bod(inp1,inp2(,inp3))
%                       inp1 - Zhler der Z-transformierten Strecke
%                       inp2 - Nenner der Z-transformierten Strecke
%          Fhrende Nullen in "inp1" und nachgestellte Nullen in "inp2" werden
%          ausgewertet! Nenner(1)=1 wird erwartet.
%       b, neue MATLAB Syntax: LTI-Modell
%          [out1(,out2)] = bod(inp1(,inp2))
%                       inp1 - LTI-Modell der Z-transformierten SISO Strecke
%   ==> Entsprechend der Syntax erfolgt die Ergebnisausgabe in zwei Versionen:
%       a, alte Syntax: out1 - Reglerzhlerpolynom (fallende Z-Potenzen)
%                       out2 - Reglernennerpolynom (    "         "    )
%                       out3 - Filterzhler
%                       out4 - Filternennerpolynom (    "         "    )
%       b, neue Syntax: out1 - LTI-Modell Regler
%                       out2 - LTI-Modell Filter
%   ==> (,out3,out4) bzw. (,out2) ist nur bei Stroptimierung sinnvoll.
%   ==> inpx = inp3(alte Syntax) bzw. inp2(LTI) = Abschaltung der Menfhrung:
%       inpx(1) - bel./1 = Optimierung: Fhrungsverhalten/Strverhalten
%       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 Frungs- / zunehmend bei Stroptimierung)
%       inpx(5) = Ordnung des zu >>berechnenden<< Reglerzhleranteils
%                 (Gesamtzhlerordnung: Kompensation+Berechnung)
%       inpx(6) - bel./1/2 = Unterdrckung von Ausgabetexten: nein/ja(K)/ja(W+K)
%                 (W - Warnung; K - Kommentar)
%       inpx(7) - bel./1 = Sonderfall gedmpfter D-Anteil einschlieen: nein/ja
%       inpx(8) - 0/1/2 - Auswahl aus drei Lsungen der kub. Gl. fr b1 (inpx(7)=1)
%       inpx(9) - bel./1 =Lsung ber quadratische Gleichung - Auswahl: V02/V01
%                 (Erfahrungswert: V02 als Verstrkung des offenen Kreises)
%       inpx(10)- bel./1/2 =Lsung ber Suchalgorithmus - Anfangswertvorgabe:
%                 maschinell/von Hand PID-Reglerparameter/von Hand Reglernullstelle
%       inpx(11)...inpx(13)=ohne Bedeutung   fr 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 Ergnzung fehlender Elemente von inpx: Null

% SCHLSSELWORTE: Digitales Betragsoptimum, Digitaler Regler von PID-Typ,
%		            Regelschleife, Stroptimierung, Fhrungsoptimierung,
%		            Fhrungsgrenfilter, Symmetrisches Optimum,
%		            Menfhrung wahlweise, Polkompensation wahlweise, LTI

%  ***** Voraussetzungen: Vorhandensein von bodg, vobo ****
%                     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 fr elektrische Antriebe.
%                     VDE-Verlag GmbH Berlin und Offenbach 1996
%                     ISBN 3-8007-1847-2
%-------------------------End German----------------------------------------------

FT01='Input argument error';
FT02='Old syntax: input of the coefficients of the plant as vectors';
FT03='Denominator(1)=1 expected!';
FT04='Input argument control_vector(inpx) - invalid dimension';
FT05='Input argument control_vector(inpx) - element 4 or/and 5 incorrect!';
FT06='Pole compensation impossible: that plant features no poles';
FT07='Calculation impossible with that setting!';
FT08='Not helpful - two degrees of freedom necessary for delayed input signals';
FT09='Incorrect input for a choice of a solution regarding b1';
FT10='Choice of a solution for b1: complex controller parameters are incorrect';
FT11='Optimization for delayed inputs impossible - too less plant structure';
FT12='Delayed input signals (PID): not enough poles for automatic initial values';
FT13='Delayed input signals: controller order > PID not prepared!';
WT01='Warning: order of the input vectors > 9 in fact?';
WT02='Warning: use of the algorithm for systems with feed through!';
WT03='Warning: dead time > 3 in fact?';
WT04='Warning: this setting enables no optimal solution!';
WT05='Warning: number of poles < number of poles intended to be compensated!';
WT06='Warning: automatic increase of the controller order due to complex poles!';
WT07='Warning: make sure integral part before mixing point of the delayed signal!';
WT08='Warning: with damped derivative part - useful in special cases only!!';
WT09='Warning: such high controller order > 4 useful?';
WT10='Warning: only  -1 < b1 < 0  corresponds with a delay time constant';
KT01='The plant features following poles:';
KT02='Compensation of the pole: ';
KT03='Compensation of the poles: ';
KT04='The cubical equation regarding b1 results in the following solutions:';
KT05='Delayed input signals (PI): solutions of the open loops gain: (V01 V02)';
KT06='Delayed input signals (PID): search algorithm needs initial values!';
KT07='Initial value of the gain of the open loop: (V01 V02)';
KT08='Initial values of the search algorithm (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-Ergnzung
if STEUV(4)<0|rem(STEUV(4),1)~=0 error(FT05); end
if STEUV(5)<0|rem(STEUV(5),1)~=0 error(FT05); end
%fr restliche STEUV-Prfung kann Eingabeprfung 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('Optimization: undelayed(any)/delayed(1) =');
else ans=STEUV(1); end
if isempty(ans)|ans~=1 ans=0; end
opa(1)=ans;
if   MUEF~=0 ans=input('Pole compensation: no(any)/yes(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 %Stroptimierung
if sden==1 error(FT06); end
pz=roots(d);
if opa(1)==0 pz=sort(pz); end %STEUV: Unterschied str- o. frungsopt.
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('no(any)/yes(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('No'); else disp('Yes'); 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('Controller integral part: yes(any.)/no(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 mglich
if   MUEF~=0 ans=input('Order of the numerator part to be computed (Enter=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 %2 Freiheitsgrade fr Stropt. ntig
if((ans==0)&(opa(1)==0)&(DKo>1))
if MUEF~=0 ans=input('"damped" derivative part: no(any)/yes(1) =');
else ans=STEUV(7); end
if isempty(ans)|ans~=1 ans=0;
else if STEUV(6)~=2 disp(WT08); end; end %Nur in Sonderfllen 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             %Fhrungsoptimierung
if opa(5)==0          %ohne Dmpfung 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 Dmpfung 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('Choice of a solution (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 unzulssig
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   % Stroptimierung
if D==1 %stroptimaler 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 %Lsungen offener Kreis
if MUEF~=0 ans=input('Choice: V02(any)/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 %stroptimaler PID-Regler
if STEUV(6)~=2&STEUV(6)~=1 disp(KT06); end
if MUEF~=0 ans=input('Choice: automatic(any)/manual(1) =');
else ans=STEUV(10); if ans==2 ans=1; end; end %Zusammenfassung fr MUEF=0
if isempty(ans)|ans~=1 ans=0; end
if ans==1
if MUEF~=0
ans=input('Choice: PID-Setting(any)/ControllerZero(1) =');
else ans=STEUV(10)-1; end
if isempty(ans)|ans~=1 ans=0; end
if ans==1
if MUEF~=0 ans=input('Controller zero = '); 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('Initial value VR= '); else ans=STEUV(11); end
anfv(1)=ans*siz(1,1);
if MUEF~=0 ans=input('Initial value d1= '); else ans=STEUV(12); end
anfv(2)=ans*anfv(1);
if MUEF~=0 ans=input('Initial value 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('Choice: V02(any)/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);                                    %Reglerzhler
if opa(1)==1 [j,i1]=size(erg);                          %Stroptimierung
out3(i1)=0; out3(1)=rzs/out1(1);           %Filterzhler
out4=erg/out1(1);                          %Filternenner
else out3=1; out4=1; end                                %Fhrungsoptimierung
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
```