function [out1,out2,out3,out4]=eez_gen(inp1,inp2,inp3)
% eez_gen calculation of discontinuous controllers via Finite Time Settling
%         method
%Differences between eez_gen.m and eez.m:
% - controller parameters: exclusively for delayed inputs, assuming: step response
%                          no zero compensation, one integral part;
%                          for desired zero compensation please use eez.m
% - suboptimal design:     choice of 3 possible approaches, related to y_max
% - check conditions:      applied for controller and pre-filter
% - inclusiv of:           optional digital behaviour illustration based on dstep
%Input/Output versions automatocally detected and same as eez.m:
% possibility a: numerator and denominator
%                inp1 - numerator   z-transformed plant (leading zeros detected)
%                inp2 - denominator z-transformed plant (leading "1" expected)
%                out1, out3 - numerator   controller and pre-filter
%                out2, out4 - denominator controller and pre-filter
% possibility b: LTI transfer function
%                inp1 - LTI transfer function z-transformed plant
%                out1 - LTI transfer function controller
%                out2 - LTI transfer function pre-filter
%Automatical menu switch off via optional definition of control vector CV
%                            = inp 3 (input version a,) or inp2 (input version b,)
%    CV(1): {1,2,3}       Cases;   assusuming y(n)=ymax; 1: n=0 2: n=1; 3: n=2
%    CV(2): {0,1,2,3,(4)} Adel;    suboptimal design (x>0) for delayed inputs
%                                  all cases; y(n)=...=y(n+x); x=4 for case 1 only
%    CV(3): {0,1,2,3,4}   Aundel;  suboptimal design (x>0) for undelayed inputs
%                                  case 1 only, i.e. y(0)=...=y(0+x);
%    CV(4): {0,1}, messages for check of controller and pre-filter; 0-off, 1-on
%    CV(5): {0,1}, figures for controlled (x) and actuating (y) variable both
%                  delayed (disturbance) and undelayed (set point) step response;
%                  0-off, 1-on
%    CV - if uncompletely automatically includes completion via "zero-elements"
%    CV: menu driven inputInstead of definition of control vector ; "any" == 0

% KEY WORDS: Finite Time Settling control, discontinuous control, control loop,
%            optimization regarding delayed input signals, optimization regarding
%            delayed input signals, pre-filter design, optionally menu driven, 
%            optionally LTI syntax, suboptimal, equalizing controller outputs

%  Copyright (c) 2014. All Rights Reserved.
%                     Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%                     https://eeiema.et.tu-dresden.de/ae2_files/ae_0e.htm
FT00='Abort: wrong case definition detected!';
FT01='Abort: case1 wrong A definition detected!';
FT02='Abort: case2 wrong A definition detected!';
FT03='Abort: case3 wrong A definition detected!';
FT04='Input arguments: general error';
FT05='Input argument CV (control_vector): incorrect dimension';
FT06='Wrong case selection';
FT07='Wrong number of suboptimal delayed steps';
FT08='Wrong number of suboptimal undelayed steps';
FT09='Wrong message outout switch';
FT10='Wrong figure outout switch';

if (nargin<1)||(nargin>3)||(nargin==2&&isa(inp2,'lti'))...
             ||(nargin==3&&isa(inp3,'lti')),   error(FT04); end
if isa(inp1,'lti')          %input is: LTI model
   if (nargin==3)||~issiso(inp1)||~isdt(inp1), error(FT04); end
   if ~isa(inp1,'tf'), inp1=tf(inp1); end;   [zz,nz,Tab]=tfdata(inp1,'v'); LTI=1;
   if nargin>1, MUEF=0; CV=inp2; else MUEF=1; CV(5)=0; end
else                        %input is: numerator + denominator polynomial
   if nargin<2, error(FT04); end;                 zz=inp1; nz=inp2;        LTI=0;
   if nargin>2, MUEF=0; CV=inp3; else MUEF=1; CV(5)=0; end
end
[znum,snum]=size(CV); if znum>1||snum>5, error(FT05); end
if snum<5, CV(5)=0; end %default completion if CV incomplete delivered
if MUEF
 Ca=input('Case (y_max) number {1,2,3] = ');
 A=input('Number A_del of suboptimal steps for delayed inputs {1,2,3,(4)} = ');
 Aun=input('Number A_und of suboptimal steps for undelayed inputs {1,2,3,4} = ');
 MS=input('Messages for controller and filter check {0,1} = ');
 FS=input('Figures for x and y delayed and undelayed step responses {0,1}= ');
else Ca=CV(1); A=CV(2); Aun=CV(3); MS=CV(4); FS=CV(5);
end
switch Ca  case {1,2,3}, otherwise error(FT06), end
switch A   case {0,1,2,3},   case {4}, if Ca~=1, error(FT07), end, 
                             otherwise error(FT07), end
switch Aun case {0,1,2,3,4}, otherwise error(FT08), end
switch MS  case {0,1}, otherwise error(FT09), end
switch FS  case {0,1}, otherwise error(FT10), end

[zer,pol,gai]=tf2zpk(zz,nz);
N=length(zer); k=length(zz)-1-N; M=length(find(pol));
if A==0, Ca=1; %Cases Ca only differ for suboptimal design, i.e. A>0; 
               %without this setting 'Ca=1' wrong results would arise
 if MS disp('Adel=0: no suboptimal delayed steps --> Cases devoid of influence!');
 end; end
D=M+A; B=N+k-1+A; Var=M+N+k+A;
%Var = number of necessary equations 
%    = number of matrix rows = number of matrix columns
%    = number of unknown to be determined
%    = (D + B + 1) - A,                 %"+1" due to VR
if length(zz)-k<Var, Nvec=zz(1+k:end); Nvec(Var)=0;
else                 Nvec=zz(1+k:Var+k); end; Nvec=Nvec/zz(k+1); %div by gain
dnz=conv([1 -1],nz);        %nz: potentially zero-filling at the end does'n mind
if length(dnz)<Var+1, dMvec=dnz; dMvec(Var+1)=0;          %"+" due to m0=1, 
else                  dMvec=dnz(1:Var+1); end;            %but VeRi start is: -m10
MaLe=tril(toeplitz(Nvec));                                %left part
if k>1, MaLe(k:end,:)=MaLe(1:end-(k-1),:);  MaLe(1:(k-1),:)=0; end%shift necessary
MaRi=tril(toeplitz(dMvec(1:end-1)));                      %right part
Mat=[MaLe(:,1:D+1) MaRi(:,1:B)];
VeRi=-dMvec(2:end)';
%selection of applicable multiplicators: consideration of A and Ca
Nmax=Nvec; if length(Nvec)<6, Nmax(6)=0; end %in general n5 necessary
n1=Nmax(2); n2=Nmax(3); n3=Nmax(4); n4=Nmax(5); n5=Nmax(6);
VM1=0;VM2=0;VM3=0;VM4=0;VM5=0;   MM2=0;MM3=0;MM4=0;MM5=0;   SM3=0;SM4=0;SM5=0;
if A~=0
   switch Ca
      case 1  %applicable for y(0)=y_max
          VM1=-n1; VM2=-n2+n1^2; VM3=-n3+2*n1*n2-n1^3;
                                 VM4=-n4+2*n1*n3-3*n1^2*n2+n2^2+n1^4;
          if     A==1, VM4=0; VM3=0; VM2=0;
          elseif A==2, VM4=0; VM3=0;
          elseif A==3, VM4=0;
          elseif A>4,  error(FT01); end
      case 2  %applicable for y(1)=y_max
                       VM2=-n2;  VM3=-(n3-n1*n2); VM4=-(n4-n1*n3+n1^2*n2-n2^2);
                       MM2=-n1;  MM3=-(n2-n1^2);  MM4=-(n3-2*n1*n2+n1^3);
          if     A==1, VM4=0;    VM3=0;    MM4=0;  MM3=0;
          elseif A==2, VM4=0;              MM4=0;
          elseif A>3, error(FT02); end
      case 3  %applicable for y(2)=y_max
                       VM3=-n3; VM4=-(n4-n1*n3); VM5=-(n5-n1*n4+n1^2*n3-n2*n3);
                       MM3=-n2; MM4=-(n3-n1*n2); MM5=-(n4-n1*n3+n1^2*n2-n2^2);
                       SM3=-n1; SM4=-(n2-n1^2);  SM5=-(n3-2*n1*n2+n1^3);
          if     A==1, VM5=0;   VM4=0;           MM5=0; MM4=0;  SM5=0; SM4=0;
          elseif A==2, VM5=0;                    MM5=0;         SM5=0;
          elseif A>3,  error(FT03); end
      otherwise error(FT00)
   end
end
% Disp_Ca_A_3MultGr=[Ca A VM1 VM2 VM3 VM4 VM5 MM2 MM3 MM4 MM5 SM3 SM4 SM5],
zN(1+5+(k-1):Var+5+(k-1),1)=Nvec; %leading zeros, easy solution avoiding [] problem
KV1=zN(5:Var+4); KV2=zN(4:Var+3); KV3=zN(3:Var+2); KV4=zN(2:Var+1); KV5=zN(1:Var);

%delete A coloumns as parameter related will be set, starting at column "1+Case"
for i=1:A, Mat(:,1+Ca)=[]; end                     % automatically column move up
LV1=Mat(:,1); LV2=Mat(:,2);                        % save column vectors 1 and 2
%"correct" left side matrix columns; multiplicators not applicable set zero before
LV1 =LV1 +KV1*VM1+KV2*VM2+KV3*VM3+KV4*VM4+KV5*VM5; % action for cases 1, 2 and 3
LV2 =LV2         +KV2*MM2+KV3*MM3+KV4*MM4+KV5*MM5; % action for cases    2 and 3
Mat(:,1)=LV1; Mat(:,2)=LV2;                        % set column vectors 1 and 2
if Var>2       % include EXCEPTION A=0 AND Var==2, i.e. M=1, N=0, k=1 --> D=1, B=0
  LV3=Mat(:,3);                                      % save column vector 3
  LV3 =LV3               +KV3*SM3+KV4*SM4+KV5*SM5;   % "correct" action for case 3
  Mat(:,3)=LV3;                                      % set column vector 3
end
Par=Mat\VeRi;   %shall be faster+more accurate as: inv(Mat)*VeRi; %compute results
%result analysis
VR=Par(1)/zz(1+k);                    %necessary consideration of dead time
b_voI=ones(1,B+1); d_vec=ones(1,D+1); %general preparation
switch Ca
  case 1, dx=[VM1 VM2 VM3 VM4];                            %predetermined formulas
          if A~=0, d_vec(1,2  :A+1)=dx(1:A); end           %set values
                   d_vec(1,2+A:D+1)=Par(2    :1+D-A)/Par(1);%computed values
                   b_voI(1,2  :B+1)=Par(2+D-A:end);
  case 2, d1=Par(2)/Par(1);
          dx=[d1*MM2+VM2  d1*MM3+VM3  d1*MM4+VM4];         %predetermined formulas
          d_vec(1,2)=d1;                                   %computed value
          d_vec(1,3:2+A)=dx(1:A);                          %set values, A>0 locked
          d_vec(3+A:1+D)=Par(    3:1+D-A)/Par(1);          %more computed values
          b_voI(1,2:B+1)=Par(2+D-A:end);
          
  case 3, d1=Par(2)/Par(1);    d2=Par(3)/Par(1);
          dx=[d2*SM3+d1*MM3+VM3 d2*SM4+d1*MM4+VM4 d2*SM5+d1*MM5+VM5];%predet.form.
          d_vec(2)=d1;         d_vec(3)=d2;                %computed values
          d_vec(1,4:3+A)=dx(1:A);                          %set values, A>0 locked
          d_vec(4+A:1+D)=Par(    4:1+D-A)/Par(1);          %more computed values
          b_voI(1,2:B+1)=Par(2+D-A:end);
end
ZRN=VR*d_vec; b_vec=conv([1 -1],b_voI); %configuration: one vector each
if MS %output controller check?
   ZP=conv(ZRN,zz); NP=conv(b_vec,nz); ZP_L=length(ZP); NP_L=length(NP);
   if ZP_L>NP_L, NP(ZP_L)=0; elseif ZP_L<NP_L, ZP(NP_L)=0; end % same dimension
   ProSu=num2str(sum(abs(ZP+NP)),'%E');
   disp(['Controller design result check: ZR*ZS+NR*NS = ' ProSu]),
end
%pre-filter design
nz_A=nz;
if length(nz)<5, nz_A(5)=0; end
m1=nz_A(2); m2=nz_A(3); m3=nz_A(4); m4=nz_A(5); %settings independent of Aun
A_vec=[1   -m1   -m2+m1^2   -m3+2*m1*m2-m1^3   -m4+2*m1*m3-3*m1^2*m2+m2^2+m1^4];
zfil_vec=A_vec(1:Aun+1);
VFi=1/(Par(1)*sum(Nvec)*sum(zfil_vec));  %VR*Vz=Par(1) - see above and VR=VR_del
ZFN=VFi*zfil_vec;  nfil_vec=d_vec;       %configuration: one vector each
if MS %output filter check?
   LimFi=num2str(sum(ZFN)/sum(nfil_vec),'%E');
   disp(['Filter design result check: Limes GFi(z-->1) = ' LimFi]),
end
if FS %output 4 step responses?
  % x: controlled variable;  y: actuating variable
  % check x and y for delayed inputs
  LN1=length(ZP); N1=1; N1(LN1)=0; figure(1), grid on, dstep(ZP,N1),  %==ZR_del*ZS
  h1=gcf(1); set(h1, 'Name','G_y_del = y / d');
  ZNPs=conv(zz,b_vec);
  LN2=length(ZNPs); N2=1; N2(LN2)=0; figure(2),grid on,dstep(-ZNPs,N2)%==ZS*NR_del
  h2=gcf(2); set(h2, 'Name','G_x_del = x / d');
  % check x and y for undelayed inputs
  VR_undel=1/(zz(1+k)*sum(Nvec)*sum(zfil_vec));
  Gg_F_del=VR_undel*conv(zfil_vec,zz);
  LN3=length(Gg_F_del); N3=1; N3(LN3)=0;
  figure(3), grid on, dstep(Gg_F_del,N3)                    %==VR_undel*A_undel*ZS
  h3=gcf(3); set(h3, 'Name','G_x_undel = x / w');
  Gg_y_und=VR_undel*conv(zfil_vec,nz);
  LN4=length(Gg_y_und); N4=1; N4(LN4)=0;
  figure(4), grid on, dstep(Gg_y_und,N4)                    %==VR_undel*A_undel*NS
  h4=gcf(4); set(h4, 'Name','G_y_undel = y / w');
end
if   LTI,  out1=tf(ZRN,b_vec,Tab);  out2=tf(ZFN,nfil_vec,Tab);
else       out1=ZRN; out2=b_vec;    out3=ZFN; out4=nfil_vec;   end
