function [Reg]=bod_es(name,T,Par,SV,AWV,OPT,P1,P2)
%BOD_ES  Digital Magnitude Optimization for control structure defined in "struc..."
%        [Reg]=bod_es(name,T,Par(,SV)(,AWV)(,OPT)(,P1,P2))
%        name=file name as string; file has to define structure of plant and
%             control (see form file: struc999.m and example: struc000.m)
%        T=   scalar definig sampling time;
%             column vector in case of parameter variants
%        Par= parameter row vector - numerical values of plant parameters,
%             sequence and meaning has to be defined in file "name";
%             matrix in case of parameter variation;
%        Reg= controller parameter - row vector, sequence and meaning has to
%                                                be defined in file "name";
%             matrix in case of parameter variation;
%             if SV(2)=1: plants zeros and poles in z-domain
%
%	 (,SV=control vector for operation mode of bod_es.m - optionally!
%	   SV(1)=0 internal MATLAB gradient estimation (default),
%              please note: in this case OPTIONS(9)=1 illegal!
%           =1 external gradient estimation for strucxxx based on bod_esp.
%              valid only for products of parameters of power "one"
%      SV(2)=0 optimization based on fsolve (default),
%           =1 no optimization, but determination of zeros and poles
%              of z-transformed plant (without additional dead time)
%              column_1=zero/pole(0/inf); column_2=respective number;
%              starting with column_3 numerical values)
%    (,AWV   = external initial values - optionally!; default: use of rand.m)
%	 (,OPT   = OPTIONS vector, see: optimization toolbox manual - optionally!)
%    (,P1,P2 = transfer of any parameters to STRUCxxx - optionally!)
%
% HINTS: *  preconditions: toolbox optimization  as well as bod_esg.m
%           trans.m, minreal.m, (and bod_esp.m but only if SV(1)=1)
%        *  In case of several parameter variants length(T) and number
%           rows inside argument Par have to be conform

% KEY WORDS: Digital Magnitude Optimum, state control structure,
%            cascade control structure, single control loop,
%            delayed input signals, undelayed input signals,
%            set point filter, parameter variation,
%            any dead time, optional mean value measurement,
%            symbolical plant definition (definition in STRUC000.m),
%            symbolical control structure definition (def. in STRUC000.m)

% Copyright (c) 2018-2023 Dr. Gert-Helge Geitner, TU Dresden,
%                     https://wwwpub.zih.tu-dresden.de/~geitner/ae_1_0e.htm
% Digital Magnitude Optimum reference: Geitner, G.-H.:
%                    Entwurf digitaler Regler für elektrische Antriebe.
%                    VDE-Verlag GmbH Berlin und Offenbach 1996
%                    ISBN 3-8007-1847-2

narginchk(3,8);     % obsolete versions: error(nargchk(3,8,nargin));
if nargin==3 SV=[]; AWV=[]; OPT=[]; P1=[]; P2=[];
elseif nargin==4 AWV=[]; OPT=[]; P1=[]; P2=[];
elseif nargin==5 OPT=[]; P1=[]; P2=[];
elseif nargin==6 P1=[]; P2=[];
elseif nargin==7 P2=[]; end
if isempty(SV) SV=[0 0]; end; if length(SV)==1 SV(2)=0; end
if isempty(OPT) OPT=optimset('TolFun',10^-30,'TolX',10^-20); end
%if SV(1)==0 GRAD=[]; else GRAD='bod_esp'; end %old call
i=length(T); [BerS,An]=size(Par); Reg=[];
if BerS~=i error('Unequal number of rows for T and Par!'), end

for iBerS=1:BerS  % BerS=computation loops

% load user defined plant structure:
[ZVe,NVe,Q,IN,OU,MT,NV]=feval(name,0,0,P1,P2);
if iBerS==1
   disp('Assumed plant parameter sequence in Par :'),
   disp(NV),      % non-recurring user info
end

% unpack name vector of plant parameters and assign numerical values:
 i=0; n0=0; k=length(NV);
 while i<k n0=n0+1; i=i+1; j=1; vtext=0;  %increase "i" by "2" because of: "; "
       while NV(i)~=';' vtext(j)=NV(i); i=i+1; j=j+1; end
       eval([vtext '=Par(' num2str(iBerS) ',' num2str(n0) ');']);
       i=i+1; end
if n0~=An
  error('Wrong plant parameter name vector NV or matrix Par!'),
end

% unpack numerator vectors of particular blocks each:
 s=1; i=0; n1=0; k=length(ZVe);
 while i<k n1=n1+1; i=i+1; j=1; vtext=0;
       while ZVe(i)~=';'  vtext(j)=ZVe(i); i=i+1;  j=j+1; end
       eval(['sz' num2str(n1) '=' vtext ';']); end

% unpack denominator vectors of particular blocks each:
 i=0; n=0; k=length(NVe);
 while i<k n=n+1; i=i+1; j=1; vtext=0;
       while NVe(i)~=';'  vtext(j)=NVe(i); i=i+1;  j=j+1; end
       eval(['sn' num2str(n) '=' vtext ';']); end
if n~=n1 error('Unequal number of plant blocks numerators / denominators!'), end

% generate unconnected diagonal system:
a=[];b=[];c=[];d=[];
for i=1:n
    [at,bt,ct,dt]=tf2ss(eval(['sz' num2str(i)]),eval(['sn' num2str(i)]));
    [a,b,c,d]=append(a,b,c,d,at,bt,ct,dt); end
% generate block connections (Q=connection matrix; IN=inputs; OU=outputs):
[i1,j1]=size(Q); [i2,j2]=size(IN); [i3,j3]=size(OU); [i4,j4]=size(MT);
if i1~=n|i2~=1|j2~=1|i3~=1
 error('Wrong block connection matrix, input or output vector'), end
if i4~=j3|j4~=2
 error('Wrong info about partial plants transformation!'), end
[ac,bc,cc,dc]=connect(a,b,c,d,Q,IN,OU); % connected system
[ap,bp,cp,dp]=minreal(ac,bc,cc,dc);     % ===> state elimination
[ZL,NL]=ss2tf(ap,bp,cp,dp,1);           %transfer functions for OU/IN
for i=1:j3                              %it has to be valid: ii=j3, i=1
    [zl,nl]=minreal(ZL(i,:),NL);
    eval(['[zz' num2str(i) ',nz' num2str(i) ']=trans(zl,nl,T(iBerS),' ...
            num2str(MT(i,1)) ',' num2str(MT(i,2)) ');']);
    if SV(2)==1  %display of zeros / poles in z-domain instead optimization
                 %delete zeros at the end before computation of roots:
           eval(['zz=zz' num2str(i) ';']); j=length(zz); ro=0;
           while zz(j)==0; j=j-1; end; for jj=1:j ro(jj)=zz(jj); end
           ro=roots(ro); j=length(ro);
           HiM(2*i-1,1)=0; HiM(2*i-1,2)=j;  %code: zeros, number
           if j~=0 HiM(2*i-1,3:3+j-1)=ro(1:j,1)'; end
           eval(['nz=nz' num2str(i) ';']); j=length(nz); ro=0;
           while nz(j)==0; j=j-1; end; for jj=1:j ro(jj)=nz(jj); end
           ro=roots(ro); j=length(ro);
           HiM(2*i,1)=inf; HiM(2*i,2)=j;    %code: poles, number
           if j~=0 HiM(2*i,3:3+j-1)=ro(1:j,1)'; end
    % zz1,nz1 % if enabled then display of z-transfer functions
    end
zl=0; nl=0; end

if SV(2)==1
     % display of plants zeros and poles in z-domain only:
   if iBerS==1
      disp('Zeros(column1=0) and poles(column1=inf) of defined');
      disp('plant transfer functions in z-domain (column2=number): ');end
   disp(HiM);            % cycle control via result display
   Reg=[Reg;HiM]; HiM=0;

else % further computation including optimization

%subsuming of single transformed transfer functions via state space description
%because of optional non-whole-number dead time and variants of measurement and
%coupling for each single transfer function:

% generate unconnected system - here variable initial values:
    a=[];b=[];c=[];d=[];at=[];bt=[];ct=[];dt=[];ac=[];bc=[];cc=[];dc=[];Q=0;
% provide an additional block as system input for function connect.m:
    eval(['zz' num2str(j3+1) '=1;']); eval(['nz' num2str(j3+1) '=1;'])
% unconnected system generation - here state description + append.m:
    for i=1:j3+1
    [at,bt,ct,dt]=tf2ss(eval(['zz' num2str(i)]),eval(['nz' num2str(i)]));
    [a,b,c,d]=append(a,b,c,d,at,bt,ct,dt); end
% connected system - here connection matrix + generation of vectors IN/OU:
    for i=1:j3+1 Q(i,1)=i; Q(i,2)= j3+1; end
    Q(j3+1,2)=0; IN=[j3+1]; OU=[1:j3];
[ac,bc,cc,dc]=connect(a,b,c,d,Q,IN,OU); % connected system
[a,b,c,d]=minreal(ac,bc,cc,dc);         %===> unconditional state elimination
if j3~=1 b(1,j3)=0; d(j3,j3)=0; end     % adaptation to connect if nblock > 1

% load control structure
[ZVe,NVe,QQ,BME,BMA,BM,NV,BRF]=feval(name,1,OU,P1,P2);
% check of user defined matrices:
[i1,j1]=size(BME);[i2,j2]=size(BMA);[i4,j4]=size(BM);[i5,j5]=size(QQ);
if i1~=i2|i2~=i4-1|j1~=j2|j2~=j4
      error('Wrong matrix BME and/or BMA and/or BM !'), end; TV=BM(i4,:);
if sum(TV)~=0 error('Wrong last row in Matrix BM !'), end
if BRF<1|BRF>j1|rem(BRF,1)~=0 error('Wrong BRF value!'), end
if max(max(BMA))>j3
   error('BMA output block number in z-domain > length(OU) !'), end
if max(max(BME))>i5
   error('BME input block number in z-domain > number of rows in QQ !'), end
j=0; for i=1:length(NV) if NV(i)==' ' j=j+1; end, end %j=number of unknown
if j~=sum(sum(BM))/j4
   error('Wrong vector NV or matrix BM !'), end

PA=sum(BM);      % determine number of parameters per variant (column)
PA_BRF=PA(BRF);  % number of parameters for selected variant
awv=length(AWV); % length of external initial value vector
if awv~=0&PA_BRF~=awv
                 error('Wrong number of external initial values!'), end
[i,j]=min(BM);   % i: minima (0!) per column; j: related row number
SZ=j(BRF)-1;     % fsolve call loop counter
enr=1;           % result number
erg=[];          % initila value result vector

TOL=sqrt(eps);   %default: unlike use of function kaskade.m
for iS=1:SZ
   TAZ=BM(iS,BRF); % transient number
   % pack fsolve arguments - now optional / obsolete version unconditionally
   knr=length(erg);
   EZW=[BRF PA_BRF TAZ iS knr];             % pack single values
   if   awv==0 AW=rand(1,TAZ); %*10,        % variable number of initial values
   else AW=0; AW(1:TAZ)=AWV(knr+1:knr+TAZ), % internal / external
   end
   for i=1:knr EZW(2,i)=erg(i); end
   BMEA=[BME BMA];
%q=fsolve('bod_esg',AW',OPT,GRAD,a,b,c,d,EZW,QQ,BMEA,ZVe,NVe,NV);%obsolete call
q=fsolve('bod_esg',AW',OPT,a,b,c,d,EZW,QQ,BMEA,ZVe,NVe,NV,TOL,SV); 
   for j=1:TAZ erg(enr)=q(j); enr=enr+1; end
       if iBerS==1&iS==1
       disp('Interpretation of result columns:'),
       disp(NV),
       end
   disp(erg);             % cycle control via result display
end

Reg(iBerS,:)=erg(1,:);
end % end of optimization

end % end of iBerS loop
