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

simu.m
```function [AM,t,vx,vy,vw2,vz,vq,vxm,vqm]=simu(fz,fn,rz,rn,zs1,sdef,wv,zv)
%SIMU     step response of a single control loop - simulation and quality criteria
%         [AM,t,vx,vy,vw2,vz,vq,vxm,vqm]=simu(fz,fn,rz,rn,zs1,sdef,wv,zv)
%Outputs:        AM - analysis matrix=[hun t(hun) tar hue t(hue) taus;
%                                     RF1 RF2    RF3 RF4 RF5    0    ]
%  hun,hue,tar,taus - undershoot resp. overshoot; rise time resp. transient time
%  t(hun),t(hue)    - point in time of the occurrence of undershoot resp. overshoot
%  RF1,RF2          - quadratic resp. absolute linear sampling integral of error
%  RF3              - enhanced quadratic sampling integral of error
%                     (additional weighted square of the controller output)
%  RF4,RF5          - ITAE (Integral of Time And Error) -sampling integral of error
%                     resp. IT2AE (quadratic time weighting)
%  t,vx,vy          - time, actual value, controller output
%  vw2,vz,vq        - filtered setpoint, disturbance variable, uncontrolled plant
%  vxm,vqm          - x resp. q mean values(MW), if mean value measuring
%Inputs:
%  fz,fn,rez,rn     - filter(f) and controller(r): numerator, denominator: z-domain
%  zs1              - polynomial between controller and disturbance action z-domain
%  sdef             - model of the continuous plant - optionally result of trans.m
%  wv,zv  - setpoint resp. disturbance variable vector, length=sampling_step_number
%Graphic representation: default colours depend on MATLAB version
%  cyan - setpoint(w2);     yellow - controller output; green - uncontrolled plant
%  blue - setpoint(20% wv); red - disturbance variable; magenta - actual value
%  also valid for mean values

% KEY WORDS: simulation, quality criteria, single control loop

%                     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

%-----------------------------German----------------------------------------------
%SIMU     Simulation und Gte der berggsfkt. fr einschleifigen Regelkreis
%         [AM,t,vx,vy,vw2,vz,vq,vxm,vqm]=simu(fz,fn,rz,rn,zs1,sdef,wv,zv)
%Ausgnge:       AM - Auswertematrix=[hun t(hun) tar hue t(hue) taus;
%                                     RF1 RF2    RF3 RF4 RF5    0    ]
%  hun,hue,tar,taus - Unter- bzw. berschwingweite; An- bzw. Ausregelzeit
%  t(hun),t(hue)    - Zeitpunkt des Auftretens von hun bzw. hue
%  RF1,RF2          - Quadratische- bzw. Betragslineare Abtast-Regelflche
%  RF4,RF5          - ITAE bzw. IT2AE (Zeitquadrat) - Abtast-Regelflchen
%  t,vx,vy          - Zeit, Istwert, Stellgre
%  vw2,vz,vq        - gefilterter Sollwert, Strgre, ungeregelte Strecke
%  vxm,vqm          - x- bzw. q-Mittelwerte(MW); wenn Mittelwertmessung
%Eingnge:
%  fz,fn,rez,rn     - Filterzhler,-nenner; Reglerzhler,-nenner im Z-Bereich
%  zs1              - Zhlerpol. zw. Regler u. Strgreneingr.   " Z-  "
%  sdef             - Def. der kont. Strecke (Ausgang von "trans.m")
%  wv,zv            - Sollwert- bzw.Stoergrenvektor (Lnge=Abtastschrittzahl)
%Graphik(Farben von MATLAB Version abhngig):
%  cyan - Sollwert(w2);     gelb - Stellgre; magenta - Istwert(auch MW);
%  blau - Sollwert(20% wv); rot  - Stoergre; grn    - ungeregelte Strecke(auch MW)

% SCHLSSELWORTE: Simulation, Gte, Einschleifiger Regelkreis

%                     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
%--------------------End German---------------------------------------------------

[i1,Df]  = size(fz);  [i2,Bf] =size(fn);
[i3,D]   = size(rz);  [i4,B]  =size(rn);
[i5,i6]  = size(sdef); N2=sdef(1,1); M=sdef(2,1); MWM =sdef(3,4);
TA=sdef(3,1); m=sdef(3,2); kode=sdef(3,3);
[i7,ios1]= size(zs1); [i8,md] =size(wv);  [i9,j]=size(zv);
if ((i1+i2+i3+i4+i7+i8+i9)~=7)|(i5~=3)|...
sum(abs(sdef(1,2+sdef(1,1):end)))~=0|sum(abs(sdef(2,2+sdef(2,1):end)))~=0
error('Incorrect inputs'); end
if (md~=j) error('Input signal vectors of different length !'); end
for i=1:N2 zs2(i)=sdef(1,1+i); end; for i=1:M ns2(i)=sdef(2,1+i); end
iof=max(Df-1,Bf-1); ior=max(D-1,B-1); ios1=ios1-1;
ans=input('Output_step_size/sampling_time: standard(10:1)   yes(any)/no(1): ');
if isempty(ans)|ans~=1 Tk=10;
else
Tk=input('Input Output_step_size/sampling_time= ');
if ((Tk<1)|(rem(Tk,1)~=0)) error('Incorrect input'); end
end
gb=input('Graphic representation: standard=100% scaling   yes(any)/no(1): ');
nov=ones(1,4);
if isempty(gb)|gb~=1 gb=0;
else j=1;
Te='Input weighting factors (x,y,w2,z) curve0 = ';
while j<5
Te(41)=48+j; Te=setstr(Te); i=input(Te);
if      ((i<0)|(i>1)) disp('Incorrect input');
elseif  i<0.1         disp('Not useful');
elseif  isempty(i)    j=5;
else nov(j)=i; j=j+1; end; end
end
if MWM
MWAZ=input('Graphic representation of mean values x and q:  no(any)/yes(1): ');
if isempty(MWAZ)|MWAZ~=1 MWAZ=0; end
else MWAZ=0; end
T=TA/Tk;j=fix((m*10^6)/(T*10^6)); m=rem(m*10^6,T*10^6)/10^6;
% Die Erweiterung um "10^6" soll richtige Ergebnisse im wahr-
% scheinlichen Zahlenbereich sichern, da "fix" und "rem" fr
% Zahlen < "1" getrennt betrachtet Fehler liefern knnen:
% z.B.: fix(0.7/0.1)=6 und rem(0.7,0.1)=0.1
mdk=fix(j/Tk); Tkk=rem(j,Tk);  %Totzeiten: bez. Abtastzeit bzw. kont. Zeitauflsung
[z0,zn]=trans(zs2,ns2,T,m,kode-MWM);   N2=length(z0); % l=ceil(m/T), m=l*T-m
M =length(zn); %M: totzeitunabhngig
if MWM==1
[z0m,znm]=trans(zs2,ns2,TA,m,kode); N3=length(z0m);%im Normalfall N3=N2+1; aber:
M3=length(znm); end%trans: Ordnungsgleichheit
j=0; while abs(z0(j+1))<(1e-14) j=j+1; end            %MWM: Totzeit j bleibt gleich
for i=(1+j):N2 zz(i-j)=z0(i); end               %fr filter.m um Totzeit verschieben
if MWM for i=(1+j):N3 zzm(i-j)=z0m(i); end; end
ios2=max(N2-j,M)-1; if MWM ios2m=max(N3-j,M3)-1; end  %Vektorlnge Anfangswerte
if iof ~=0 zif(iof)=0;   end
if ior ~=0 zir(ior)=0;   end
if ios1~=0 zis1(ios1)=0; end
if ios2~=0 zis2(ios2)=0;
ziur(ios2)=0; end; x=0;
if MWM&ios2m~=0 zis2m(ios2m)=0;
ziurm(ios2m)=0; end;
y2v(mdk+1)=0; sev(mdk+1)=0; vxv(Tkk+1)=0; vqv(Tkk+1)=0;
vx(1)=0;  vy(1)=0;  vw2(1)=0; vz(1)=0; vw(1)=0; vxw(1)=0; vq(1)=0;
if MWM vxm(1)=0; vqm(1)=0; xm=0; qm=0; end
for i=1:md
if   iof~=0  [w2,zff] =filter(fz,fn,wv(i),zif);  zif=zff;
else         [w2]     =filter(fz,fn,wv(i));      end
if MWM  xw=w2-xm;
else    xw=w2-x; end
if ior~=0  [y,zfr]  =filter(rz,rn,xw,zir);       zir=zfr;
else  [y]      =filter(rz,rn,xw);           end
if ios1~=0 [y1,zfs1]=filter(zs1,1,y,zis1);       zis1=zfs1;
else  [y1]     =filter(zs1,1,y);            end
y2=y1-zv(i);    y2v(1)=y2; y2=y2v(1+mdk); %einschreiben,auslesen,schieben
se=wv(i)-zv(i); sev(1)=se; se=sev(1+mdk);
for j=1:mdk  y2v(mdk+2-j)=y2v(mdk+1-j); %Totzeit bezglich Abtastzeit
sev(mdk+2-j)=sev(mdk+1-j); end
for j=1:Tk
if ios2~=0 [x,zfs2]=filter(zz,zn,y2,zis2); zis2=zfs2;
[q,zfur]=filter(zz,zn,se,ziur); ziur=zfur;
else  [x]     =filter(zz,zn,y2);
[q]     =filter(zz,zn,se);      end
vxv(1)=x;     x=vxv(1+Tkk); %einschreiben,auslesen,schieben
vqv(1)=q;     q=vqv(1+Tkk);
for ii=1:Tkk  vxv(Tkk+2-ii)=vxv(Tkk+1-ii);     %Totzeit bezglich konti-
vqv(Tkk+2-ii)=vqv(Tkk+1-ii); end %nuierlicher Zeitauflsung
vx(2+j+(i-1)*(Tk+1)) =x;
vy(1+j+(i-1)*(Tk+1)) =y;
vw2(1+j+(i-1)*(Tk+1))=w2;
vz(1+j+(i-1)*(Tk+1)) =zv(i);
vw(1+j+(i-1)*(Tk+1)) =wv(i);
vxw(1+j+(i-1)*(Tk+1))=xw;
vq(2+j+(i-1)*(Tk+1)) =q;
if MWM vxm(1+j+(i-1)*(Tk+1))=xm;
vqm(1+j+(i-1)*(Tk+1))=qm; end
end
if MWM
if ios2m~=0 [xm,zfs2m]=filter(zzm,znm,y2,zis2m); zis2m=zfs2m;
[qm,zfurm]=filter(zzm,znm,se,ziurm); ziurm=zfurm;
else        [xm]      =filter(zzm,znm,y2);
[qm]      =filter(zzm,znm,se);      end
end
vx(3+j+(i-1)*(Tk+1)) = vx(2+j+(i-1)*(Tk+1));% Treppenfunktn ermglichen
vy(2+j+(i-1)*(Tk+1)) = vy(1+j+(i-1)*(Tk+1));%       "            "
vw2(2+j+(i-1)*(Tk+1))=vw2(1+j+(i-1)*(Tk+1));%       "            "
vz(2+j+(i-1)*(Tk+1)) = vz(1+j+(i-1)*(Tk+1));%       "            "
vw(2+j+(i-1)*(Tk+1)) = vw(1+j+(i-1)*(Tk+1));%       "            "
vxw(2+j+(i-1)*(Tk+1))=vxw(1+j+(i-1)*(Tk+1));%       "            "
vq(3+j+(i-1)*(Tk+1)) = vq(2+j+(i-1)*(Tk+1));%       "            "
if MWM
vxm(2+j+(i-1)*(Tk+1))=vxm(1+j+(i-1)*(Tk+1));%    "            "
vqm(2+j+(i-1)*(Tk+1))=vqm(1+j+(i-1)*(Tk+1));end% "            "
end
vw2(3+j+(i-1)*(Tk+1))=vw2(2+j+(i-1)*(Tk+1));% Zeitvektor auffllen
vz(3+j+(i-1)*(Tk+1)) = vz(2+j+(i-1)*(Tk+1));%     "          "
vw(3+j+(i-1)*(Tk+1)) = vw(2+j+(i-1)*(Tk+1));%     "          "
if MWM  xw=w2-xm;
else    xw=w2-x;  % w2-Brechnung nicht ntig, mit Filter keine Auswertung
end
if ior~=0  [y,zfr]  =filter(rz,rn,xw,zir);
else    [y]      =filter(rz,rn,xw); end
vxw(3+j+(i-1)*(Tk+1))=xw;
vy(3+j+(i-1)*(Tk+1))=y;
if MWM vxm(3+j+(i-1)*(Tk+1))=xm;
vqm(3+j+(i-1)*(Tk+1))=qm; end
szma=max([max(wv) max(zv)]); szmi=min([min(wv) min(zv)]);
if szma<0 szma=0; end
if szmi>0 szmi=0; end
wabs=max([szma abs(szmi)]); vwn = vw/wabs*100;
vw2n=vw2/wabs*100;
vxn = vx/wabs*100;
vzn = vz/wabs*100;
i=2; while sdef(1,i)==0; i=i+1; end; Vp=sdef(1,i);
vqn = vq/(wabs*Vp)*100;
if MWM vqmn=vqm/(wabs*Vp)*100;
vxmn=vxm/wabs*100; end
snmi=szmi/wabs*100; snma=szma/wabs*100;
wema=max([max(vw2n) max(vxn)]); wemi=min([min(vw2n) min(vxn)]);
if (snma-wema)<0 yma=snma-(snma-wema); else yma=snma; end
if (snmi-wemi)>0 ymi=snmi-(snmi-wemi); else ymi=snmi; end

stma=max(vy); if stma<0 stma=0; end
stmi=min(vy); if stmi>0 stmi=0; end
wabs=max([stma abs(stmi)]); vyn = vy/wabs*100;
stnma=stma/wabs*100; stnmi=stmi/wabs*100;
if stnma>yma yma=stnma; end
if stnmi<ymi ymi=stnmi; end
if gb vxn=vxn*nov(1); vyn=vyn*nov(2); vw2n=vw2n*nov(3); vzn=vzn*nov(4);
if MWM vxmn=vxmn*nov(1); end; end
axis([0 md ymi*1.05 yma*1.05])
t=-1:Tk-1; t(1)=0;           % Darstellung von Treppenkurven ermglichen
for i=1:md-1 for j=1:Tk+1 t((Tk+1)*i+j)=t(j)+Tk*i; end, end
t(1+md*(Tk+1))=t(md*(Tk+1))+1; t(2+md*(Tk+1))=t(1+md*(Tk+1)); t=t/Tk;
%4.2c:
%plot(t,vxn,'c11',t,vyn,'c5',t,vw2n,'c12',t,vzn,'c2',t,vwn*0.2,'c10',...
%     t,vqn,'c3');
%hher:
if ~MWAZ plot(t,vqn,'g',t,vyn,'y',t,vw2n,'c',t,vzn,'r',t,vwn*0.2,'b',...
t,vxn,'m');
else     plot(t,vqn,'g',t,vyn,'y',t,vw2n,'c',t,vzn,'r',t,vwn*0.2,'b',...
t,vxn,'m',t,vxmn,'m',t,vqmn,'g');
end
xlabel('time/samplig_time'); ylabel('per_cent'); grid;
disp('Stop: please press any key.'); pause
vxR=vx; vyR=vy;           %Retten wegen nachfolgender Abfrage zum Ausgangssetzen
if MWM vxmR=vxm; end
ans=input('Outputs = simulation_data(SD) or drawings(DR) ? SD(any)/DR(1): ');
if isempty(ans) ans=0; end
if ans==1 vx=vxn; vy=vyn; vw2=vw2n; vz=vzn; vq=vqn;
if MWM vxm=vxmn; vqm=vqmn; end; end
ans=input('Analyses of the response behaviour no(any)/yes(1): ');
if isempty(ans) ans=0; end
if ans~=1 AM=0; if ~MWM vxm=0; vqm=0; end
else TvP=1;  %Position des letzten Zeitwertes vor dem Sprung
while ((wv(TvP)==0)&(zv(TvP)==0)&(TvP<md)) TvP=TvP+1; end; bte=0;
if zv(TvP)~=0 bfu=0; prwe=zv(TvP);
if sum(abs(wv))~=0|sum(abs(zv))~=abs(sum(zv)); bte=1; end
else          bfu=1; prwe=wv(TvP);
if sum(abs(zv))~=0|sum(abs(wv))~=abs(sum(wv)); bte=1; end
end
Tv=t(1+(Tk+1)*(TvP-1)); %normierter Zeitwert bei Sprung
if bte disp([07 07 'Prepared for a single stepped response only !']);
else
vxn=vxn/nov(1); vyn=vyn/nov(2); vw2n=vw2n/nov(3); vzn=vzn/nov(4);
ans=input('Range of tolerance: standard=1%   yes(any)/no(1): ');
if isempty(ans)|ans~=1 toba=1;
else i=0;
while i==0 toba=input('Input range of tolerance (in %) = ');
if  toba>10|toba<0|isempty(toba) disp('Illegal value');
else i=1; end; end
end;
disp('Choice of a sampling integral of error - (0=all,Enter=no): ');
if ((brf<0)|(brf>5)|(rem(brf,1)~=0)) error('Incorrect input'); end
if isempty(brf) brf=6; end
if ((brf==5)|(brf==0))
ans=input('Weighted controller output: standard=0.1   yes(any)/no(1): ');
if isempty(ans)|ans~=1 lampda=0.1;
else lampda=input('Input - controller output weighting = ');
if lampda<0|lampda>1|isempty(lampda) error('Incorrect input'); end;
end; end
if MWM
MWAW=input('Quality criteria by instantaneous_value(any)/mean_value(1): ');
if isempty(MWAW)|MWAW~=1 MWAW=0; else vxn=vxmn/nov(1); end
else MWAW=0; end
disp('Quality criteria of the selected calculation intreval result in:');
if prwe<0 vxn=-vxn; vwn=-vwn; vzn=-vzn; end
[hun,t1P]=min(vxn); t1=t(t1P)-Tv;
if abs(hun)<0.01 hun=0; t1=0; disp('Undershoot represents 0 %');
else disp('Undershoot in unit % at time_t/sampling_time:');
disp([hun t1]); end
[hue,t2P]=max(vxn); t2=t(t2P)-Tv;
if bfu hue=hue-100; vgv=vwn; pep=1+toba/100; mep=1-toba/100;
else                vgv=vzn; pep=toba/100;   mep=-pep; end
if hue<0.01 hue=0; t2=0; disp('Overshoot represents 0 %');
else disp('Overshoot  in unit % at time_t/sampling_time:');
disp([hue t2]); end
i=t1P; an=2+md*(Tk+1); j=an+1; %Strsprung: eintauchen statt verlassen erkennen
while i<an
if  (vxn(i+1)<=(vgv(i+1)*mep)) i=i+1; else j=i+1; i=an; end; end
if j==an+1 tar=inf; else tar=t(j)-Tv; end
disp('Rise_time/sampling_time = '); disp(tar);
i=0; j=an;
while i<an
if ((vxn(an-i)<=(vgv(an-i)*pep))&(vxn(an-i)>=(vgv(an-i)*mep)))
i=i+1; else j=an-i+1; i=an; end; end
if    j==an+1; taus=inf; elseif j==an taus=0; else taus=t(j)-Tv; end
disp('Transient_time/sampling_time = '); disp(taus);
AM(1)=hun; AM(2)=t1; AM(3)=tar; AM(4)=hue; AM(5)=t2; AM(6)=taus;
if brf~=6
RF1=0; RF2=0; RF3=0; RF4=0; RF5=0; vyn=vyR/prwe;
if MWAW vxw=(vw-vxmR)/prwe; else vxw=(vw-vxR)/prwe; end
tn=t; k=0;                                                  %Doppelung streichen
for i=1:Tk+1:an tn(i-k)=[]; vxw(i-k)=[]; vyn(i-k)=[]; k=k+1; end;
Anz=an-k;                  %verringerte Werteanzahl wegen gestrichener Doppelung
Bg=1+(TvP-1)*Tk; %Korrektur fr Berechnungsbeginn bei verzgertem Eingangssignal
if (brf==2)|(brf==0)
for i=Bg:Anz RF2=RF2+abs(vxw(i)); end; end
if (brf==1)|(brf==0)
for i=Bg:Anz RF1=RF1+(vxw(i))^2; end; end
if (brf==4)|(brf==0)
for i=Bg:Anz RF4=RF4+(tn(i)-Tv)*abs(vxw(i)); end; RF4=RF4*TA; end %Echtzeitbezug
if (brf==5)|(brf==0)
for i=Bg:Anz RF5=RF5+(tn(i)-Tv)^2*abs(vxw(i)); end; RF5=RF5*TA^2; end %   "
if (brf==3)|(brf==0)
for i=Bg:Anz RF3=RF3+((vxw(i))^2)+lampda*((vyn(i)-vyn(Anz))^2); end; end;
vRF=[RF1 RF2 RF3 RF4 RF5]*(TA/Tk); %konstante lineare Delta-t-Auflsung beachten
if brf==0 disp('Integrals of error represent:');
format short e; disp(vRF); format; brf=5;
else disp('The selected integral of error represents:'); disp(vRF(brf)); end
for i=1:brf AM(2,i)=vRF(i); end
end
end
end
```