function [RZ,FZ]=bod_gen(SZ,CV,AW,options)
%bod_gen Discontinuous Amplitude Optimum for any single control loop
%        [RZ (FZ)]=bod_gen(SZ,CV) - with control vector CV
%        [RZ (FZ)]=bod_gen(SZ)    - OR optionally menu-driven
%        [RZ,(FZ)]=bod_gen(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 Amplitude Optimum (BOD) to calculate the
%coefficients of a discontinuous controller RZ (and optionally pre-filter FZ
%if applicable); without any internal pole compensation - if desired: previous
%external realization necessary!
% 
% SZ:	                   packed plant in z-domain (LTI system)
% CV=[D B (LL) (IP) (PF)]: control vector otherwise menu-driven monitoring
% D:	                   order controller numerator
% B:	                   order controller denominator
% (LL):  default=0    switch: controller type [general(0)/LeadLag(1)]
% (IP):  default=1    switch: first order integral part [0=off / 1=on]
% (PF):  default=0    switch: pre-filter [0=off / 1=on]; order==D
% RZ:                 packed controller (LTI system) in z-domain
% (FZ):               packed pre-filter (LTI system) in z-domain (delayed inputs)

% KEY WORDS: Discontinuous Amplitude Optimum, BOD, single control loop,
%            optimization, undelayed / delayed input signals, general controller,
%            delayed derivative part, lead and leg controller.

%  ***** preconditions: existence of function bod_gen_eq.m ****
%  Copyright (c) 2018. All Rights Reserved.
%                     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_0.htm
%  Reference regarding Discontinuous Amplitude 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 Lead/Lag controller: Roberto Zanasi, Riccardo Morselli:
%                     Disrete inversion formulas for the design of lead and
%                     lag discrete compensators. proceedings ECC 2009, pp.5069/74
%                     ISBN 978-963-311-369-1; F.T. Brown: Engineering System
%                     Dynamics 2. edition, pp.551/552

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 5 elements!';
FT04='Plant SZ in z-domain has to be a LTI system!';
FT05='D, B, LL, IP and PF have to be of type scalar!';
FT06='LL, IP and PF have to be of type "switch", i.e. "0" or "1"!';
FT07='D and B have to be of type natural number!';
FT08='Integral part or pre-filter not prepared for phase lead / lag controller!';
FT09='Only first (2 dof) or second order (3 dof) lead / lag controller prepared!';

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=');
     LL=input('switch: controller type [general(0)/LeadLag(1) = ');
     IP=input('switch: first order integral part [0=off / 1=on] = ');
     PF=input('switch: pre-filter (order==D) [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>5 error(FT03); end
     D=CV(1); B=CV(2); LL=0; IP=1;     PF=0;
     if     col==3 LL=CV(3);
     elseif col==4 LL=CV(3); IP=CV(4);
     elseif col==5 LL=CV(3); IP=CV(4); PF=CV(5); end
end
%common input argument check
if ~isa(SZ,'lti') error(FT04); end
if ~isscalar(D)||~isscalar(B)||~isscalar(LL)||~isscalar(IP)||~isscalar(PF)...
   error(FT05); end
if ~(LL==0||LL==1)||~(IP==0||IP==1)||~(PF==0||PF==1)...
   error(FT06); end
if D<0||rem(D,1)~=0||B<0||rem(B,1)~=0 error(FT07); end

if LL==1 %phase lead/lag controller
   if IP==1||PF==1 error(FT08); end             %neither prefilter nor integrator
   if ~((D==1&&B==1)||(D==2&&B==2)) error(FT09); end %only 1. or 2. order
end

if nargin<3
   if ~LL
      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
   else  beta=rand(1);                        %initial values
         if D==2 AW(3,1)=(1+2*beta)/3; end   %LL 2. order controller (3 dof)
         AW(2,1)=beta; AW(1,1)=1+2*beta;      %LL 1. order controller (2 dof)
   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_gen_eq(x,SZ,D,B,IP,PF,LL),AW,options);

% *** Attention: x return: column vector ***
[zz,nz,T]=tfdata(SZ,'v');     %concerning sampling time T;'v':otherwise cell array
if ~LL
   rz=x(1:(D+1))';               %D order controller numerator;  D=0: gain only
   rn=x((D+1):(D+1+B))';         %B order controller denominator;B=0: minimum form
   % ***
   rn(1)=1;                          %b0=1; Integral part extra!
   if IP rn=conv(rn,[1 -1]); end     %with integral part
   RZ=tf(rz,rn,T,'Variable','z^-1'); %packed controller in z-domain
   if  PF                            %optimization for delayed inputs
       fn=rz/rz(1);                  %division by gain
       fz=sum(fn);                   %final value theorem
       FZ=tf(fz,fn,T,'Variable','z^-1'); %packed controller in z-domain
   else
       FZ=tf(1,1,T,'Variable','z^-1');   %optimization for undelayed inputs
   end
else
   alpha=x(1); beta=x(2);
   if D==1  rz=[alpha  -(alpha-1)]; %first order Zanasi: identically with
            rn=[beta   -(beta-1) ]; %first order approximation of differentiation
   else     gamma=x(3);             %second order: three unknown quantities (dof=3)
            rz=[alpha  -(alpha-1+gamma) gamma]; %first+second order diff.
            rn=[beta   -(beta-1+gamma+(beta-alpha)/3)  (gamma+(beta-alpha)/3)];
   end
   RZ=tf(rz,rn,T,'Variable','z^-1'); %packed controller in z-domain
   FZ=tf(1,1,T,'Variable','z^-1');   %dummy as output arg.
end
