function [REG,FIL,A2,FLAG] = bo_d_pltgen(regtype,PLT,Ttot,app,X0,optim)
%function [REG,FIL(,A2,FLAG)]=bo_d_pltgen(regtype,PLT(,Ttot,app,X0,optim))
%Optimization of continuous controllers: Magnitude Optimum - delayed inputs
% ---> this function applies LTI definitions of transfer functions
% Plant:        PLT=     LTI plant transfer function inclusive of:
%                        gain, numerator polynomial, denominator polynomial
%                        BUT exclusive of dead time and its approximation
% Controller:   regtype={2,3,4}={PI,PID,PID2}
%               REG=    LTI controller transfer function
%                       GC(p)=[1+p*(TN)+p^2*(TV*TN)+p^3*(TQ*TN)]/[p*(TR0)],
%                       since definition VR=TN/TR0 is possible and p is
%                       Laplace variable, an alternative reconfiguration 
%                       would be possible:GC(p)=VR*[1+1/(p*TN)+p*TV+p^2*TQ]
% Pre-filter:   FIL=    LTI pre-filter transfer function
%                       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"
% optim:   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
[nump,denop]=tfdata(PLT,'v');
Vs=nump(end); nump=nump/Vs;   %calculation normalized to Vs is advantageous
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
%Check for dead time:
if     nargin<3  NO_Tt=1;            %no dead time
elseif nargin==3 
       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>3 
       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>6 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

XX=[];
if exist('X0','var')&&~isempty(X0) AW=X0;  % [TN,TR0,TV,TQ]
   else AW=100*rand(1,regtype); end, %initial values
if exist('optim','var') OPT=optim; else OPT=[];   end
  %StepToleranz=TolX;   FiniteDifferenceStepSize=FinDiffRelStep
  %if regtype==2 %equation system may be too easy for 'trust-region-dogleg'
  %FSdefault=optimset('fsolve');
  %OPT=optimset(FSdefault,'TolFun',1e-40,'FinDiffRelStep',sqrt(eps),...
  %                 'TolX',1.0e-06,...
  %                 'Algorithm','levenberg-marquardt');
  %end
warning OFF BACKTRACE  % see "help warning": control the display of the
         %file and line number that produced a warning when it is displayed
[XX,A2,FLAG]=fsolve(@(XX) bo_d_pltgen_eq(XX,regtype,nump,denop),AW,OPT);
%XX, % XX = [TN   TR0  TV  TQ  q5]
      
TN=XX(1);  TR0=XX(2)*Vs;             % Algorithm is based on T0 rated by Vs
if     regtype==2 REG=tf([TN  1],[TR0  0]);
                  FIL=tf(1,[TN  1]);
elseif regtype==3 TV=XX(3);
                  REG=tf([TV  TN  1],[TR0  0]);
                  FIL=tf(1,[TV  TN  1]);
elseif regtype==4 TV=XX(3); TQ=XX(4);
                  REG=tf([TQ  TV  TN  1],[TR0  0]);
                  FIL=tf(1,[TQ  TV  TN  1]);
else   disp('Input regtype~={2,3,4} improper'); end
XXL=length(XX);
if XXL==5 q5=XX(5); if q5~=0  
       if    q5==-1 warning('Basic condition B0=A0=1 violated'),
       elseif q5==1 warning('At least one parameter zero for condition 1'),
       elseif q5==2 warning('At least one parameter zero for condition 2'),
       elseif q5==3 warning('At least one parameter zero for condition 3'),
       elseif q5==4 warning('At least one parameter zero for condition 4'),
       else         warning('Improper error code q5');
end, end, end
warning ON BACKTRACE     
end

function q=bo_d_pltgen_eq(XX,regtype,nump,denop)
TN=XX(1);  TR0=XX(2);  rede=[TR0  0];  % XX = [TN   TR0  TV  TQ]
if     regtype==2 renu=[TN  1]; 
elseif regtype==3 TV=XX(3);  renu=[TV  TN  1];
elseif regtype==4 TV=XX(3);  TQ=XX(4);  renu=[TQ  TV  TN  1];
else   disp('Input regtype~={2,3,4} improper'); end
B= flip(nump);                            % Gg_num==Zs
A1=flip(conv(renu,nump));
A2=flip(conv(rede,denop));                % Gg_den==Nr*Ns+Zr*Zs (Nfi=Zr)
MLA=max([length(A1)  length(A2)]);
A1(MLA+1)=0; A2(MLA+1)=0;
A=A1+A2;                                  % Gg_den
ML=max([length(B)  length(A)]);
B(ML+1)=0; A(ML+1)=0;
if ML<9 B(10)=0; A(10)=0; end % so limited till condition 4 ONLY!

%condition 0
q0=B(1)-A(1); if q0~=0 q5=-1; end  % error code condition 0

% condition 1 (P-, I-, D-) --- must be calculated in any case
if ML<3 q5=1, end           % error code condition 1
q(1)= ( B(2)^2-2*B(3)*B(1) )...
     -( A(2)^2-2*A(3)*A(1) );

if regtype==4          % condition 4 (PID2)
if ML<9 q5=4, end           % error code condition 4
q(4)=  ( B(5)^2-2*B(4)*B(6)+2*B(3)*B(7)-2*B(2)*B(8)+2*B(1)*B(9) )...
     - ( A(5)^2-2*A(4)*A(6)+2*A(3)*A(7)-2*A(2)*A(8)+2*A(1)*A(9) );
end

if regtype<5&&regtype>2  % condition 3 (PID)
if ML<7&&regtype>1 q5=3, end  % error code condition 3
q(3)=  ( B(4)^2-2*B(3)*B(5)+2*B(2)*B(6)-2*B(1)*B(7) )...
     - ( A(4)^2-2*A(3)*A(5)+2*A(2)*A(6)-2*A(1)*A(7) );
end

if regtype<5&&regtype>1  % condition 2 (PI-, PD-)
if ML<5&&regtype>1 q5=2; end  % error code condition 2
q(2)=  ( B(3)^2-2*B(2)*B(4)+2*B(1)*B(5) )...
     - ( A(3)^2-2*A(2)*A(4)+2*A(1)*A(5) );

end
if exist('q5','var') q(5)=q5; end % possible feedback, would cause:
% Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square
% systems; using Levenberg-Marquardt algorithm instead.
end