frfn.m

MFD－多变量系统频域设计工具

function [numd,fit,k]=frfn(w,f,charl,cdend,weight)
%FRFN   Constructs a Reversed-Frame Normalizing controller.
%       [NUMD,fit,K]=FRFN(W,F,CHARL,CDEND,WEIGHT) calculates
%       a compensator numerator matrix, NUMD, given:
%       -  A system to be compensated ,F, in MVFR form,
%       -  An associated frequency vector, W,
%       -  A set of desired frequency responses for the characteristic gain
%              loci of the compensated system. These are specified
%              in the columns of the matrix CHARL. One column per output.
%       -  A common denominator transfer function matrix, CDEND,
%       -  An optional argument, WEIGHT, a matrix of column vectors
%             of the frequency response of the weightings.
%             One column per output of the system.
% If fsize(F) is [m,n] then fsize(WEIGHT) and fsize(CHARL) are both [1,m]
%                    m must be less than or equal to n
%       FRFN(W,F,CHARL,CDEND,WEIGHT) returns a numerator
%       matrix, NUMD, for the pre-compensator. CDEND^(-1)*NUMD.
%
%       [NUM,fit]=FRFN(W,F,CHARL,CDEND,WEIGHT) also returns the
%       relative error, fit, as an indicator of how well the compensator
%       meets the desired response,
%
%       [NUM,fit,K]=FRFN(W,F,CHARL,CDEND,WEIGHT) also returns the
%       resulting compensator evaluated over W as the MVFR matrix, K.

%       J-M. Boyle 2 Sep 1987. Revised by J.M.Maciejowski, 3 Jan 1988.
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd
%       Modified by M.P. Ford 5th October 1987
%       to accept a non-common denominator Cdend

nargs = nargin; nargo = nargout;
error(nargchk(4,5,nargs));

[mf, nf] = fsize(w, f);       % Find the size of the component matrices.
if mf>nf
   error('System must have at least as many inputs as outputs.')
end
lw = length(w);

% Define 'i' and establish the number of frequency points.
jay = sqrt(-1);

%  Calculate the SVD of G(s).
[yf, sf, uf] = fsvd(w, f);
yfct  = ftrn(w, yf);

disp('FRFN  Setting-up matrices')

% calculate the matrix Q.
[cm,cn]=fsize(w,charl);
if (cm~= 1) | (cn~=mf)
   error('Fsize(CHARL) needs to be 1 by m, where fsize(F) is [m,n]')
end
charl=fdiag(w,charl);
q = fmulf(w, fmulf(w, yf, charl), yfct);   % Q(s) = Y(s).charl(s).Y(s)

% Specify the weighting matrix.
  if nargs == 4
    weight = eyef(w,mf);
  else
    [wm,wn]=fsize(w,weight);
    if (wm~=1) | (wn~= mf)
       error('Fsize(WEIGHT) needs to be 1 by m, where fsize(F) is [m,n]')
    end
    weight = weight.*conj(weight);
    weight = fdiag(w,weight);
  end
   % weighting matrix is abs(weight)^2

% Set-up the denominator matrix.
[md, nd] = size(cdend);       % Find the order of the denominator.

% Calculate the size of the numerator matrix.
if rem(nd,nf)~=0
   error(['Columns of Denominator must be an integer',...
	  ' multiple of the system inputs'])
end
dimd=nd/nf;
numd=zeros(dimd,1);
numd(dimd)=1;
% Calculate the MVFR matrix for S(s)/D(s)
if md~=nf
   error('Denominator does not have number of rows equal to number of inputs.')
end
d = mv2fr(cdend,numd, w);

% Calculate the size of the frequency part of coefficient matrix.
a = zeros(nf*lw, dimd*mf);

% Calculate the values of the coefficient matrix at each frequency.
for m = 1:lw                   % For each frequency value
  for i=1:nf                       % Repeat for each row
    k = 0;                         % Counter for the 'order' of element.
    for j=(dimd*i):-1:(dimd*(i-1)+1)
       a(i+nf*(m-1), j) = (jay*w(m)).^k; %  All other terms are jw(i) .
       k = k + 1;
    end           % of s terms
  end             % of row
end               % of frequency points

d=finv(w,d);
				  %       -1
a = fmulf(w, d, a);               %  D(jw) .S(jw)
if nargo >= 3                   % Store S to calculate compensator
   s_save = a;
end
% Calculate the Matrix A.         %            -1
a = fmulf(w, f, a);               % A = G(jw).D (jw).S(jw)

% ^      T    -1     T
% n = ( A W A )  .( A W q )
% 'n' contains the values of the numerator coefficients.
disp('FRFN  Calculating the compensator numerator matrix')
numd = ( real( a'* fmulf( w, weight, a)))\(real(a' * fmulf(w, weight, q)));
if nargo >= 2
   n_save = numd;
end

% Sort out the rows and columns.
n2 = [];nout=[];
for j=1:nf
   for i=1:nf
      n2 = [n2;numd( (j-1)*dimd+1:j*dimd,i)];
   end
   nout = [nout;n2']; n2 = [];
end
numd=nout;

% Calculate the indicator of how good fitting is.
disp('FRFN  Calculating fit')
if nargo >= 2
   fit = norm( faddf( w, q, (-1*fmul( w, a, n_save))), 2)./(norm( q, 2));
end

% calculate the compensator.
if nargo == 3                     % if compensator required
				    %         -1
   k = fmul(w,s_save,n_save);       % K = D (s) * S(s) * N
end