facc.m

控制系统计算机辅助设计——MATLAB语言与应用(源代码）

function [ka,kb,kc,kd] = facc(w,f,index,snum,sden)
%FACC   Synthesise an approximate commutative controller.
%
%       [ka,kb,kc,kd] = facc(w,f,index,snum,sden)
%
%       f is the MVFR matrix of a square system and w is the 
%       associated frequency list. The approximate commutative 
%       controller is designed at frequency w(index). 
%       (index should be a single integer.)
%       The i'th row of snum and the i'th row of sden contain the 
%       numerator and denominator, respectively, of the i'th diagonal
%       compensator of the controller, namely the compensator applied
%       to the i'th characteristic locus.
%
%       ka,kb,kc,kd are the state space matrices of the resulting 
%       controller.

%       Dr J.M. Maciejowski,  20 November 1987. Revised 7 Dec 1987.
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd

[m,n]=fsize(w,f);
if m~=n
  error('Not a square system. (FACC)')
end

% Check consistency of dimensions:
[nr,nc]=size(snum);
if nr ~= n
  error('Incorrect number of numerators. (FACC)')
end
[dr,dc]=size(sden);
if dr~=n
  error('Incorrect number of denominators. (FACC)')
end
%ka=[]; kb=[]; kc=[]; kd=[];  % newly added statement

% First get a state space model of the diagonal part:
for i=1:n
  % Remove leading zeros:
  sni = snum(i,min(find(snum(i,:)~=0)):length(snum(i,:)));
  sdi = sden(i,min(find(sden(i,:)~=0)):length(sden(i,:)));
  % Get state space model of i'th diagonal:
  [ai,bi,ci,di]=tf2ss(sni,sdi);
  if abcpgchk(ai,bi,ci)     % Deal with pure gains
    ai=[]; bi=[]; ci=[];
  end
  ni=max(size(ai));
  ka = [ka,  zeros(max(size(ka)),ni);
	zeros(ni,max(size(ka))), ai  ];
  [br,bc]=size(kb);
  kb = [kb, zeros(br,1);
	zeros(ni,bc), bi ]; 
  [cr,cc]=size(kc);
  kc = [kc, zeros(cr,ni);
	zeros(1,cc), ci  ];  
  kd = [kd, zeros(cr,1);
	zeros(1,bc), di ];
end
% Assertion: kb,kd have n columns, kc,kd have n rows.
%
% Now find real approximations to the eigenvector frames:
[v,d]=eig(fgetf(w,f,index));
wr=align(v);
vr=align(inv(v));  % vr=inv(wr) may work better sometimes
%
% Now form series connection vr*s*wr:
kb=kb*wr;
kc=vr*kc;
kd=vr*kd*wr;
%
% Uncomment next line to ensure minimal realization:
% [ka,kb,kc,kd]=minreal(ka,kb,kc,kd);