function  Reg = bo_d_PI(regtype,Vs,nump,denop,Ttot,app)
%function Reg = bo_d_PI(regtype,Vs,nump,denop(,Ttot,app))
%Optimization of continuous controllers: Magnitude Optimum - delayed inputs
% Controller:   regtype={2}, i.e. PI only; alternative to bo_d_gen
%               via direct solution of a quadratic equation here
%               Reg=[VRa,TNa,TR0a; VRb,TNb,TR0b] - set 'a' typically useful
%                      parameter meaning: GC(p)=VR*[1+1/(p*TN)]
%                      VR=TN/TR0, p: Laplace variable, thus alternatively:
%                      GC(p)=[1+p*TN]/(p*TR0)
% Plant:        Vs=gain; nump  = numerator polynomial;
%                        denop = denominator polynomial
% Pre-filter:   pre-filter denominator = controller numerator (ex gain)
%               GFi(p)=1/[1+p*TN]
% Ttot:    optional dead-time
% app:     optional approximation code for Ttot: default 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

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,
elseif nargin>6 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:
if     regtype==2 xP=3; xD=1; % PI ==>  PS=3, DZ/DP: 1
else   error('Error: wrong input "controller type" in bo_d_PI'); 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

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);
K0=PS0^2*PS3+2*PS1*DP1;  K1=PS1*PS2-PS0*PS3;  K2=PS3+2*PS1*DZ1;
if K0==0 warndlg(['Weighting factor VR^0 not applicable, '...
                  'not a quadratic equation'],...
                 'PI special case K0'),  % unlikely to happen
    Vra=2*K1/K2;                        Vrb=Vra;
elseif K1==0 warndlg(['weighting factor VR^1 not applicable, ' ...
                  'K0 or K2 should be negative'],...
                  'PI special case K1'), % check K1
    Vra=sqrt(-K0/K2);                   Vrb=Vra;
elseif K2==0 warndlg(['Weighting factor VR^2 not applicable, ' ...
                  'not a quadratic equation'],...
                  'PI special case K2'),  % not a quadratic equation
    Vra=K0/(2*K1);                      Vrb=Vra;
else
    Vra=K1/K2*(1+sqrt(1-K0*K2/K1^2));   Vrb=K1/K2*(1-sqrt(1-K0*K2/K1^2));
end
    T0a=-2*PS1/(Vra^2+2*PS0*Vra+PS0^2); T0b=-2*PS1/(Vrb^2+2*PS0*Vrb+PS0^2);

%Definitionen für Regler mit einem I-Anteil
VRa=Vra/Vs;    VRb=Vrb/Vs;   % Vr   == XX(1)  is identical to variable z
                             % 1/T0 == XX(2)  is identical to variable y
TR0a=T0a*Vs;   TR0b=T0b*Vs;  % Algorithm is based on T0 rated by Vs 
TNa=VRa*TR0a;  TNb=VRb*TR0b; % ==Vr*T0: here rating must cancel itself out
Reg=[VRb,TNb,TR0b; VRa,TNa,TR0a];  %Reg=[VR,TN,TR0] two solution sets
   if ~d_off   %switch off disp commands - useful for nested function calls         
   disp('Result parameter sequence: [VRa, TNa, TR0a; VRb, TNb, TR0b]');
   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