function  Reg= bo_u_gen(regtype,Vs,nump,denop,Ttot,app)
%function  Reg= bo_u_gen(regtype,Vs,nump,denop(,Ttot(,app)))
%Optimization of continuous controllers: 
%Magnitude Optimum - undelayed inputs, no set-point pre-filter;
%               any consideration of a controller filter for the practical 
%               application/realization has to be "outside" as plant add-on
% plant:        Vs=gain; nump  = numerator polynomial;
%                        denop = denominator polynomial
% Controller:   regtype={1,2,3,4}={I,PI,PID,PID2}
%               Reg=    [VR,TN,T0,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)
% 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
% 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
%vorläufig nur Regler angenommen mit einem I-Anteil
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
if     nargin<5  NO_Tt=1;            %no dead time
elseif nargin==5 NO_Tt=0; app=23;    %dead time, but no approximation rule
elseif nargin==6 NO_Tt=0;
else   disp('Wrong number of input arguments'); 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:1,3,5,7; and DP:0,1,2,3
if     regtype==1 xP=1; xD=0; %I 
elseif 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_u_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
DPV=[]; for i=1:xD DPV(i)=DPX(i,nf); end     %if maxVL>=xP length=maxVL+1

Vec=[]; Mat=[]; Reg(1,5)=0;
PS0=n0; PS1=PSV(1); DP0=0.5*n0^2;
Vec=-PS1;  Mat=DP0;                                      % I ?
if regtype>1 PS2=PSV(2); PS3=PSV(3); DP1=DPV(1);         % PI ?
             Vec(2,1)=PS3; Mat(1,2)=PS0; Mat(2,1:2)=[DP1 -PS2]; end
if regtype>2 PS4=PSV(4); PS5=PSV(5); DP2=DPV(2);         % PID ?
             Vec(3,1)=-PS5; Mat(2,3)=-PS1;
             Mat(3,1:3)=[DP2 PS4 PS3]; end
if regtype>3 PS6=PSV(6); PS7=PSV(7); DP3=DPV(3);         % PID2 !
             Vec(4,1)=PS7; Mat(2,4)=-PS0; Mat(3,4)=PS2;
             Mat(4,1:4)=[DP3 -PS6 -PS5 -PS4]; end
Erg=inv(Mat)*Vec;                                     %Erg=[T0, Tn, Tp, Tq]
Reg(3)=Erg(1)*Vs;                  %T0=was rated with Vs for simplification
if regtype>1 Reg(2)=Erg(2); Reg(1)=Reg(2)/Reg(3); end %Reg=[VR,TN,T0,TV,TQ]
if regtype>2 Reg(4)=Erg(3)/Erg(2); end  %Tp was computed, TV is needed
if regtype>3 Reg(5)=Erg(4)/Erg(2); end  %Tq was computed, TQ is needed
   if ~d_off   %switch off disp commands - useful for nested function calls
disp(['Sequence of the result parameters: ' 'VR, TN, T0, TV, TQ']);
   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 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);
    DP=DP+(-1)^(x+i)*denop(i+1)*denop(2*x-i+1);
end
end