function erg=zust_bran_2(T,Tei,TM,TF,TL,SAK)
%ZUST_BRAN_2 state control structure 2: 4 coefficients; automated; here plant only,
%        erg=zust_bran_2(T,Tei,TM,TF,TL(,SAK)) - application of the Discontinuous
% Magnitude Optimum(BOD) to calculate coefficients of a discontinuous state control
% structure - example by Brandenburg: at 35(1987)7, pp. 275/83;
% quod vide: Tei=4; TM=30; TF=5; und TL=30;60;120  **** automated calculation ****
% BUT: modified control structure - now 4 parameters
%      T:	   sampling time of the discontinuous control in ms
%      Tei:	   first order time delay (model of closed current control loop) in ms
%      TM:	   motor time constant (integrator) in ms
%      TF:	   spring time constant (integrator) in ms
%      TL:	   load time constant  (integrator) in ms
%      (,SAK): optionally: switch for operation control, 
%              display of partial results / diagnostics : 0=no, nor
%              default=0;                    1=yes, display of partial results
%                                            2=yes, partial results + diagnostics
%      erg:	   results as vector [K1 K2 K3 K4]

% KEY WORDS: Discontinuous Magnitude Optimum, BOD, state control structure,
%            elastic drive train, optimization, undelayed input signals,
%            delayed input signals

%  ***** precondition: existence of function zust_bran_2gl.m ****
%  Copyright (c) 1997-2023. All Rights Reserved.
%                     Dr. Gert-Helge Geitner, TU Dresden,
%                     https://wwwpub.zih.tu-dresden.de/~geitner/ae_1_0e.htm
%  Reference regarding Discontinuous Magnitude Optimum: Geitner, G.-H.:
%                     Entwurf digitaler Regler fuer elektrische Antriebe.
%                     publisher: VDE-Verlag GmbH Berlin und Offenbach 1995
%                     ISBN 3-8007-1847-2

narginchk(5,6) % obsolete versions: error(nargchk(5,6,nargin))
if nargin==5 SAK=0; 
   else if ~(SAK==0|SAK==1|SAK==2) error('Incorrect input argument 6!'); end
end
S=[1 1 1 1];	% S=[K1 K2 K3 K4]   **** switch for weighting factors
RFV=[1 2 3 4];	% RFV=[1 2 3 4]	    **** sequence vector with [K1 K2 K3 K4]
   if sum(RFV)~=10 error('Incorrect sequence vector'), end
BM=[1 1 1; 2 3 1; 1 0 1; 0 0 1; 0 0 0];%COLUMN: number calculations per fsolve call
BRF=1;			   % **** choice of a column of matrix BM
% Definition of numerators and denominators of all blocks (transfer functions)
  n001=1;                          d001=[Tei 1];
  n002=1;                          d002=[TM 0];
  n003=1;                          d003=[TF 0];
  n004=1;                          d004=[TL 0];
% declaration of the number of blocks (transfer functions) and definition of an
% unconnected diagonal state model of the system using the above defined
% numerator and denominator polynomials:
%nblocks=4;
%blkbuild;      % because of global variable access for script files only!!!
% define an unconnected diagonal system
[a,b,c,d]=tf2ss(n001,d001);  s0n=abs('n000'); s0d=abs('d000');
for t=2:4  HUS=fix(t/100); ZEV=fix(t/10); sk=[0 HUS ZEV-HUS*10 t-ZEV*10];
           s1=char(s0n+sk); s2=char(s0d+sk);
           [at,bt,ct,dt]=tf2ss(eval(s1),eval(s2));
           [a,b,c,d]=append(a,b,c,d,at,bt,ct,dt);
end;       %  !!!Thus algorithm prepared for three ASCII positions (0...9)!!!
% definition of the connections between the blocks (transfer functions):
Q=[ 1  0  0
    2  1 -3
    3  2 -4
    4  3  0 ];    in =[1];  ou=[1 2 3 4];
[ac,bc,cc,dc]=connect(a,b,c,d,Q,in,ou); % connected system
[ap,bp,c,d]=minreal(ac,bc,cc,dc);       % state elimination
[a,b]=c2d(ap,bp,T);                     % z-transform
b(1,4)=0; d(4,4)=0;                     % rebuilding of a symmetrical
                                        % structure=number of blocks
                                        % a,b,c,d: invariable part
PA=sum(S);		% number of parameters
BA=sum(BM);		% initial values per variant (column)
if PA~=BA(BRF) error('Incorrect calculation sequence'), end
[i,j]=min(BM);		% i: Minima (0!) per column; j: associated row number
SZ=j(BRF)-1;		% loop counter for fsolve calling
enr=1;			% result number
erg=[];			% initial value
%possible OPTIONS display: optimset,
if SAK==2
OPT=optimset('Display','iter','Diagnostics','on','TolFun',10^-30,'TolX',10^-20);
%OBSOLATE OPTION Settings
%options=optimset('Display','iter','Diagnostics','on',...
%                 'TolFun',10^-30,'TolX',10^-20,'NonlEqnAlgorithm','lm');
else            %also alternative algorithm: 'gn'; BUT NOT default 'dogleg'!!!
OPT=optimset('Display','off','Diagnostics','off','TolFun',10^-30,'TolX',10^-20);
%OBSOLATE OPTION Settings
%options=optimset('Display','off','Diagnostics','off',...
%                 'TolFun',10^-30,'TolX',10^-20,'NonlEqnAlgorithm','lm');
end             %i.e. algorithm without: 'LargeScale','on'
for iS=1:SZ
 if iS~=1&SAK>0 erg, end % operation control  ****
 TAZ=BM(iS,BRF);         % transient number
 AW=ones(TAZ,1)*20;      % variable number of initial values
       %factor "20" necessary since obsolate option: 'NonlEqnAlgorithm','lm'
 %AW=randn(TAZ,1);       % variable number of initial values
 x=[];
 x=fsolve(@(x) zust_bran_2gl(x,a,b,c,d,erg,BRF,RFV,TAZ,iS),AW,OPT);
 for j=1:TAZ
     erg(enr)=x(j);
     enr=enr+1;
 end
end
