function  [Reg,A2,FLAG]= bo_d_gen(regtype,Vs,nump,denop,Ttot,app,X0,OPT)
%function [Reg(,A2,FLAG)]=bo_d_gen(regtype,Vs,nump,denop(,Ttot,app,X0,OPT))
%Optimization of continuous controllers: Magnitude Optimum - delayed inputs
% Plant:       Vs=gain; nump  = numerator polynomial;
%                       denop = denominator polynomial
% Controller:  regtype={2,3,4}={PI,PID,PID2}
%              Reg=    [VR,TN,TR0,TV,TQ]
%                      parameter meaning: GC(p)=VR*[1+1/(p*TN)+p*TV+p^2*TQ]
%                      VR=TN/TR0, p: Laplace variable, thus alternatively:
%                      GC(p)=[1+p*TN+p^2*TV*TN+p^3*TQ*TN]/(p*TR0)
% Pre-filter:  pre-filter denominator = controller numerator (ex gain)
%              GFi(p)=1/[1+p*TN+p^2*TV*TN+p^3*TQ*TN]
% Ttot:    optional dead-time
% app:     optional approximation code for Ttot: if absent app=23
%          11 ... 19 = Strejc approximation of order 01 ... 09
%          21 ... 23 = Pade         "        "   "   01 ... 03
%          31 ... 39 = Polynomial   "        "   "   01 ... 09
% X0:      optional initial vales for equation solver "fsolve"
% OPT:     optional OPTIONS for equation solver "fsolve"
% A2,FLAG: optional outputs A2 and EXITFLAG of equation solver "fsolve"
%          meaningful results: A2 as small as possible, Flag=1
% For deactivation of display commands use hidden mode "regtype+10" as type
%                                        - useful for nested function calls

%Check of inputs:
if round(regtype)~=regtype|~isscalar(regtype) 
                                        disp('Input regtype improper'); end
if regtype>10 regtype=regtype-10; d_off=1;     %realizes "disp" hidden mode
                                  else d_off=0; end
if ~isreal(Vs)||~isscalar(Vs)  disp('Input Vs improper'); end

for j=1:length(nump) if ~isreal(nump(j))||~isscalar(nump(j))
                  disp(['Input 3 nump(' num2str(j) ') improper']), end, end
for j=1:length(denop) if ~isreal(denop(j))||~isscalar(denop(j))
                 disp(['Input 4 denop(' num2str(j) ') improper']), end, end
Knz=find(denop,1,'last');            %Last vector element unequal zero
DeltaInt=length(denop)-Knz;          %plant integral terms
NumInt=DeltaInt+1;                   %controller features one integral term
%%% Previously designed for controllers with just one integral part !!!
% --------- should not be applied -------------------------------------------
%if NumInt>1             %hence, n0=1 NOT valid, but n0=...=n_(NumInt-2)==0
%    Vs=Vs/denop(Knz);       %hence, Vs contains fictive units
%    denop=denop(1:Knz);     %hence, denominator must be shortened
%    denop=denop/denop(Knz); %hence, now fictive n0==1
%end,                        n0=1; %thus, settig is general valid
% -------------------------------------------------------------------------
n0=denop(end);
if     nargin<5  NO_Tt=1;            %no dead time
elseif nargin==5 
       if isempty(Ttot) NO_Tt=1;    %[] == no dead time, no more arguments
       else NO_Tt=0; app=23; end,   %dead time, but no approximation rule
elseif nargin>5 
       if isempty(Ttot) NO_Tt=1;    %[] == no dead time, but more arguments
       else NO_Tt=0; end,      %dead time and rule, more arguments possible
elseif nargin>8 disp('Wrong number of input arguments'); %handle X0 and OPT
end, 
   if ~NO_Tt    %% plant features dead time - start
if ~isreal(Ttot)||~isscalar(Ttot)||Ttot<0  disp('Input Ttot improper'); end
if round(app)~=app|~isscalar(app)|app<0 disp('Input app improper'); end
%set chosen Ttot approximation:
if app>10&&app<20                           %approximation 1: Strejc
   ProxOrder=app-10; Tt_num=1; Tt_den=1;
   for i=1:ProxOrder Tt_den=conv(Tt_den,[Ttot/ProxOrder 1]); end,
elseif app>30&&app<40                       %approximation 3: Polynomial
   ProxOrder3=app-30; Tt_num=1; Tt_den(1,ProxOrder3+1)=1;
   for i=1:ProxOrder3 Tt_den(ProxOrder3-(i-1))=Ttot^i/factorial(i); end,
elseif app==21  Tt_num=[-Ttot  1];          %approximation 2: Pade
                Tt_den=[+Ttot  1];          %reconfiguration p^0: Koeff=1
elseif app==22  Tt_num=[(1/12)*Ttot^2  -0.5*Ttot  1];   
                Tt_den=[(1/12)*Ttot^2   0.5*Ttot  1];
elseif app==23  Tt_num=[-(1/120)*Ttot^3  0.1*Ttot^2  -0.5*Ttot  1];
                Tt_den=[+(1/120)*Ttot^3  0.1*Ttot^2   0.5*Ttot  1];
else   disp('Input app improper'); end 
%Firstly polynomial multiplication thereafter reversal of order
%Tt_num, Tt_den,  % for check ONLY
nump=conv(nump,Tt_num); denop=conv(denop,Tt_den);
   end          %% plant features dead time - end
%set polynomials in ASCENDING order now:
zf=flip(nump); nf=flip(denop);
maxVL=max(length(zf),length(nf));
zf(maxVL+1)=0; nf(maxVL+1)=0;  %Vector extension is irrelevant

% MINIMUM polynomial order for PI:  4!!! check nesssary
% regtype depended necessary orders::   PS:3,5,7;    DZ and DP:1,2,3
if     regtype==2 xP=3; xD=1; %PI
elseif regtype==3 xP=5; xD=2; %PID
elseif regtype==4 xP=7; xD=3; %PID2
else   error('Error 1 in bo_d_gen'); end
if maxVL<xP zf(xP+1)=0; nf(xP+1)=0; end      %adequate vector length
PSV=[]; for i=1:xP PSV(i)=PSX(i,zf,nf); end  %irrespective zeros and poles
DZV=[]; for i=1:xD DZV(i)=DZX(i,zf); end     %if maxVL>=xP length=maxVL+1
DPV=[]; for i=1:xD DPV(i)=DPX(i,nf); end

XX=[];
if exist('X0','var')&&~isempty(X0) AW=X0;
                              else AW=rand(1,regtype); end, %initial values
if exist('OPT','var') optim=OPT; else optim=[];   end,
  %StepToleranz=TolX;   FiniteDifferenceStepSize=FinDiffRelStep
  %if regtype==2 %equation system may be too easy for 'trust-region-dogleg'
  %FSdefault=optimset('fsolve');
  %optim=optimset(FSdefault,'TolFun',1e-40,'FinDiffRelStep',sqrt(eps),...
  %                 'TolX',1.0e-06,...
  %                 'Algorithm','levenberg-marquardt');
  %end
[XX,A2,FLAG]=fsolve(@(XX) bo_d_gen_eq(XX,n0,regtype,PSV,DZV,DPV),AW,optim);
%Definitionen für Regler mit einem I-Anteil
VR=XX(1)/Vs;          % Vr   == XX(1)  is identical to variable z
T0=1/XX(2);           % 1/T0 == XX(2)  is identical to variable y
TR0=T0*Vs;            % Algorithm is based on T0 rated by Vs 
TN=VR*TR0;            % ==Vr*T0: here rating must cancel itself out
Reg=[VR,TN,TR0];         %Reg=[VR,TN,TR0,TV,TQ]
if regtype>2 Tpr=XX(3); %Tp/T0 is identical to variable x
             Tp=Tpr*T0; TV=Tp/TN; Reg(4)=TV; end
if regtype>3 Tqr=XX(4); %Tq/T0 is identical to variable w
             Tq=Tqr*T0; TQ=Tq/TN; Reg(5)=TQ; end
   if ~d_off   %switch off disp commands - useful for nested function calls         
disp(['Sequence of the result parameters: ' 'VR, TN, TR0, TV, TQ']);
   end
end

function q=bo_d_gen_eq(XX,n0,regtype,PSV,DZV,DPV)
z=XX(1);  y=XX(2);
if regtype>2 x=XX(3); else x=0; end
if regtype>3 w=XX(4); else w=0; end
PS0=n0; PS1=PSV(1);     PS2=PSV(2); PS3=PSV(3);
DZ0=0.5; DZ1=DZV(1);    DP0=0.5*n0^2; DP1=DPV(1);
q(1)=...
   DZ0*z^2+PS0*z+DP0+PS1*y              -2*DZ0*x*y;
q(2)=...
   DZ1*z^2-PS2*z+DP1-PS3*y+DZ0*x^2-PS1*x-2*DZ1*x*y-2*DZ0*z*w        -PS0*w;
%ACHTUNG q-Anzahl bedingt berechnen
if regtype>2
PS4=PSV(4); PS5=PSV(5); DZ2=DZV(2); DP2=DPV(2);
q(3)=...
   DZ2*z^2+PS4*z+DP2+PS5*y+DZ1*x^2+PS3*x-2*DZ2*x*y-2*DZ1*z*w+DZ0*w^2+PS2*w;
end
if regtype>3
PS6=PSV(6); PS7=PSV(7); DZ3=DZV(3); DP3=DPV(3);
q(4)=...
   DZ3*z^2-PS6*z+DP3-PS7*y+DZ2*x^2-PS5*x-2*DZ3*x*y-2*DZ2*z*w+DZ1*w^2-PS4*w;
end
end 
function PS=PSX(x,nump,denop)
%assumption 1: equal length of nump and denop
%assumption 2: x and length are harmonated
PS=0;
for i=0:x
    PS=PS+(-1)^i*denop(i+1)*nump(x+1-i);
end
end
function DZ=DZX(x,nump)
%assumption: x and length of nump are harmonated
DZ=0.5*nump(x+1)^2;
for i=0:x-1
    DZ=DZ+(-1)^(i+1)*nump(x-i)*nump(x+i+2);
end
end
function DP=DPX(x,denop);
%assumption: x and length of denop are harmonated
DP=0.5*denop(x+1)^2;
for i=0:x-1
    DP=DP+(-1)^(i+1)*denop(x-i)*denop(x+i+2);
end
end