function [R, S, T]=bod_rst(SZ,CV,AW,options)
%bod_rst Discontinuous Magnitude Optimum - any single RST control loop
%        [R, S, T]=bod_rst(SZ,CV) - with control vector CV
%        [R, S, T]=bod_rst(SZ)    - OR optionally menu-driven
%        [R, S, T]=bod_rst(SZ,CV,AW,options)
%  AW      - optionally initial values, column vector (fsolve control)
%  options - optionally options, see example below    (fsolve control)
%  Please note: there isn't any verification for AW and options.
%
%Application of the Discontinuous Magnitude Optimum (BOD) to calculate the
%coefficients of discontinuous controller RST parts; without any internal 
%pole compensation - if desired: previous external realization necessary!
% 
% SZ:	                   packed plant in z-domain (LTI system)
% CV=[D B (IP)]: control vector otherwise menu-driven monitoring
% D:	                   order controller R-part
% B:	                   order controller S-part (without integrator)
% (IP):  default=1    switch: first order integral part [0=off / 1=on]
% R:                  controller R-part in z-domain
% S:                  controller S-part in z-domain
% T:                  controller T-part (gain) in z-domain

% KEY WORDS: Discontinuous Magnitude Optimum, BOD, RST-control, design

%  ***** preconditions: existence of function bod_rst_eq.m ****
%  Copyright (c) 2018-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 1996
%                     ISBN 3-8007-1847-2
%  Reference regarding RST controller: Ioan Landau, Gianluca Zito:
%                     Digital control systems - design, indentification
%                     and implementation, Springer, 2010,
%                     ISBN 978-1-84996-551-4

FT01='Control vector CV has to be a row vector but no LTI system!';
FT02='Control vector CV has to be a row vector but no matrix!';
FT03='Control vector CV should have 2 till 3 elements!';
FT04='Plant SZ in z-domain has to be a LTI system!';
FT05='D, B and IP have to be of type scalar!';
FT06='IP has to be of type "switch", i.e. "0" or "1"!';
FT07='D and B have to be of type natural number!';

narginchk(1,4) % obsolete versions:  error(nargchk(1,4,nargin))
if   nargin==1 % menu-driven input desired
     D=input('Controller order numerator D=');
     B=input('Controller order denominator B=');
     IP=input('switch: first order integral part [0=off / 1=on] = ');
else           %application of control vector CV
     if isa(CV,'lti') error(FT01); end
     [row,col]=size(CV);
     if row~=1 error(FT02); end
     if col<2||col>3 error(FT03); end
     D=CV(1); B=CV(2); IP=1;
     if       col==3   IP=CV(3); end
end
%common input argument check
if ~isa(SZ,'lti') error(FT04); end
if ~isscalar(D)||~isscalar(B)||~isscalar(IP) error(FT05); end
if ~(IP==0||IP==1) error(FT06); end
if D<0||rem(D,1)~=0||B<0||rem(B,1)~=0 error(FT07); end

if nargin<3
   AW=rand(D+B+1,1);                     %general initial values
   if     D==0&&B==0 AW(1,1)=AW(1,1)+1;  %P-  or I-controller
   elseif D==1&&B==0 AW(1,1)=AW(1,1)+1;  %PD- or PI-controller
                     AW(2,1)=-0.8*AW(1,1);
   elseif D==0&&B==1 AW(1,1)=AW(1,1)+2;  %PID-controller with ext. pole comp.
                     AW(2,1)=-0.7;       %delayed diff.-part via b1
   elseif D<4&&B==0  AW(1,1)=AW(1,1)+1;  %D=2&B=0: PID-controller
                     AW(2,1)=-1.5*AW(1,1);
                     AW(3,1)=0.7*AW(1,1);
          if D==3    AW(4,1)=-AW(4,1); end %PID2-controller
   end
end, %AW,

if nargin<4
%OPTIONS display: optimset,                     %iteration supervision
FSdefault=optimset('fsolve');
options=optimset(FSdefault,'Display','off','Diagnostics','off','TolFun',1e-18,...
                 'Algorithm','trust-region-reflective','FinDiffType','central');
%UNVALID OPTION Settings
%options=optimset('Display','off','TolFun',1e-18,'Diagnostics','off',...
%    'LargeScale','off','LevenbergMarquardt','off','LineSearchType','quadcubic',...
%    'Simplex','off','NonlEqnAlgorithm','lm','PrecondBandWidth',inf);  
%options=optimset('Display','iter','TolFun',1e-18,'Diagnostics','on',...
%    'LargeScale','off','LevenbergMarquardt','off','LineSearchType','quadcubic',...
%    'Simplex','off','NonlEqnAlgorithm','lm','PrecondBandWidth',inf);
end

x=[]; 
x=fsolve(@(x) bod_rst_eq(x,SZ,D,B,IP),AW,options);
    % FMINCON: Supplied objective function must return a scalar value
    %LB=-3*ones(D+1+B,1);
    %UB=+3*ones(D+1+B,1);
    %x=fmincon(@(x) bod_rst_eq(x,SZ,D,B,IP),AW,[],[],[],[],LB,UB,...
%x=lsqnonlin(@(x) bod_rst_eq(x,SZ,D,B,IP),AW,LB,UB,options);                                                          %[],options);

% *** Attention: x return: column vector ***
[zz,nz,T]=tfdata(SZ,'v');    %concerning sampling time T;'v': or cell array
R=x(1:(D+1))';               %D order controller R-part; D=0: gain only
  R(1)=sqrt(R(1));           % da VR positiv erzwungen: siehe ... eq_WRZL.m
S=x((D+1):(D+1+B))';         %B order controller S-part without integral
% ***
S(1)=1;                      %b0=1; Integral part extra!
if IP S=conv(S,[1 -1]); end  %with integral part if any
T=sum(R);                    %gain only, assuming open loop integral part
