sumkromat.m

剑桥大学用于线性和双线性动力学系统辨识的工具箱

function y = sumkromat(A,C,U,i)
%SUMKROMAT   Calculates the sum of block Kronecker product matrices.
%
%    SUMKROMAT(A,C,U,i) constructs the sum of a number of Kronecker products of
%      subsets of the data matrix U and an observability-type matrix for the
%      state-space system (A,B,C,D). Used in SUBID3B.
%
%   Usage:
%      Y = sumkromat(A,C,U,i)
%
%   Inputs: 
%      A := Square matrix with n rows and n columns.
%      C := Matric with p rows and n columns.
%      U := Block data matrix with 
%             U := [U(1); U(2);...; U(i); U(i+1);...; U(3i)].
%             Each U(i) has m rows and N-3i+1 columns.
%      i : = Positive integer.
%   
%   Outputs:  
%      Y := Sum of Kronecker product matrices, defined as
%                         [ I     0    ]               [ 0     0    ]
%            Y := U(i)' * [ 0  O(2i-1) ] + U(i+1)' *   [ I     0    ] + ...
%                         [ 0          ]               [ 0  O(2i-2) ]
%
%                            [ 0 0 ]  where O(i) is the observability matrix
%               + U(3i-1)' * [ 0 0 ], with index i as computed using OBMAT 
%                            [ I 0 ]  and * is the Kronecker product.
%
%             Y has (N-3i+1)2ip rows and (n+p)m columns.
%
%   See also KRON, OBMAT, SUBID3B.
%
% CUED System Identification Toolbox.
% Cambridge University Engineering Department.
% Copyright (C) 1998-2002. All Rights Reserved.

% Version 1.00, Date: 01/06/2002
% Created by H. Chen and E.C. Kerrigan.

% See (40) in Section 8.2 of N.L.C. Chui and J.M. Maciejowski. "Subspace
% identification - a Markov parameter approach". Technical Report
% CUED/F-INFENG/TR.337, December 1998, University of Cambridge, UK.

[rowA,colA] = size(A);
[rowC,colC] = size(C);
[rowU,colU] = size(U);

if colA ~= rowA,
  error('A is not a square matrix')
end
if colA ~= colC,
  error('The column numbers of A and C are not equal')
end
if (i < 1) | (i ~= fix(i)), 
  error('Input i must be a positive integer')
end

p = rowC;
n = rowA;
m = rowU/(3*i);

if m ~= fix(m)
  error('U is not divisible by 3i')
end

t = zeros(colU*2*i*p,(n+p)*m);

for j = 0:2*i-1
  
  UJ = U((i+j)*m+1:(i+j+1)*m,:)';
  
  % MOD ECK 22/6/2001
  % Description above did not agree with what is implemented in the
  % following commented region
  
  %   if j == 0,
  % 	tt = [zeros((2*i-1)*p,p+n);
  % 		 eye(p) zeros(p,n)   ];
  %   else
  % 	tt(1:(2*i-1-j)*p,:) = zeros((2*i-1-j)*p,n+p);
  % 	tt((2*i-1-j)*p+1:(2*i-j)*p,:) = [eye(p) zeros(p,n)];
  % 	tt((2*i-j)*p+1:2*i*p,:) = [zeros(j*p,p) obmat(A,C,j)];
  %   end

  if j == 2*i-1
	tt = [ zeros((2*i-1)*p,p+n);
		   eye(p)    zeros(p,n)];
  else
	tt = [ zeros(j*p,n+p);
		   eye(p)               zeros(p,n);
		   zeros((2*i-1-j)*p,p) obmat(A,C,2*i-1-j)];
  end
  
  % END MOD 
  
  t = t + kron(UJ,tt);
  
end

y = t;
 
% *** last line of sumkromat.m ***