stlsident.m

The main features of the considered identification problem are that there is no an a priori separati

 function [ sysh, info, w, xini ] = stlsident( w, m, l, opt )
% STLSIDENT - Identification by structured total least squares.
%
% [ sysh, info, wh, xini ] = stlsident( w, m, l, opt )
%
% Inputs:
% W   - given time series (#samples x #variables x #time series).
% M   - input dimension.
% L   - system lag.
% OPT - options for the optimization algorithm:
%   OPT.EXCT - vector of indices for exact variables (default []);
%   OPT.SYS0 - initial approximation: an SS system with M inputs, 
%              P:=size(W,2)-M outputs, and order N:=L*P;
%   OPT.DISP    - level of display [off | iter | notify | final];
%   OPT.MAXITER - maximum number of iterations (default 100);
%   OPT.EPSREL, OPT.EPSABS, and OPT.EPSGRAD (default 1e-4)
%     - convergence tolerances (see the help of STLS).
%
% Outputs:
% SYSH - input/state/output representation of the identified system.
% INFO - information from the optimization solver:
%   INFO.M - misfit ||W - WH||_F,
%   INFO.TIME - execution time for the STLS solver,
%   INFO.ITER - number of iterations.
%   NOTE: INFO.ITER = OPT.MAXITER indicates lack of convergence.
% WH   - optimal approximating time series.
% XINI - initial conditions under which WH are obtained. 

% Requires: STLS, CORR, DECODE_STRUCT, SYS2X, INISTATE.

% Author: Ivan Markovsky, Date: November 2004 
% Reference: I. Markovsky, J. C. Willems, S. Van Huffel, 
% B. De Moor, and R. Pintelon, "Application of structured total 
% least squares for system identification and model reduction", 
% Report 04--51, Dept. EE, K.U. Leuven, 2004

% Option names
OPTNAMES = {'disp','epsrel','epsabs','epsgrad', ...
            'maxiter','sys0','exct'};
DISPNAMES = {'off','iter','notify','final'};

% Constants
T  = size(w,1); % # of samples
nw = size(w,2); % # of variables
k  = size(w,3); % # of time series
p  = nw - m;    % # of outputs
n  = p * l;     % order of the system

% Check for parameters out of range
if nargin < 3
  error('Not enough input arguments. (W, M, and L should be given.)')
end
if (~isreal(w))
  error('Complex data not supported.')
end
if (~isreal(m) | m < 0 | m > nw | (ceil(m)-m ~=  0))
  error('Invalid number of inputs M. (0<M<size(W,2) and integer.)')
end
if (~isreal(l) | l < 0 | (ceil(l)-l ~=  0))
  error('Invalid system lag L. (0 < L and integer.)')
end

% Assign the default options
disp    = 'notify';
epsrel  = 1e-4;
epsabs  = 1e-4;
epsgrad = 1e-4;
maxiter = 100;
sys0    = [];  % to be chosen
exct    = [];  % no exact variables

% Check the user given options
if (nargin > 3)
  if isstruct(opt)
    names = fieldnames(opt);
    I = length(names);
    for i = 1:I
      ind = strmatch(lower(names{i}),OPTNAMES,'exact');
      if isempty(ind)
        warning('OPT.%s is an invalid option. Ignored.', ...
                upper(names{i}))
      else
        eval([OPTNAMES{ind} ' = opt.' names{i} ';']);
      end
    end
  else
    warning('OPT should be a structure. Ignored.')
  end
end

% Check for invalid option values
exct = unique(exct);
if any(1 > exct | nw < exct)
  error('An index for exact variable is out of range.')
end
if (length(exct)*T > m*T + n)
  error('Too many exact variables. Generically there is no solution.')
end

% Detect trivial cases
if (m == nw) 
  warning('M = size(W,2) => trivial solution SYSH = I, WH = W.');
  sysh      = ss([],[],[],eye(nw),-1);
  info.M    = 0;
  info.iter = 0;
  if nargout == 4
    xini = inistate(w,sysh);
  end
  return
end
if (T <= n) 
  warning('T < p*l => trivial solution.');
  a = diag(ones(n-1,1),-1);
  b = zeros(n,m);
  c = [w(:,m+1:end)' zeros(p,T-n)];
  d = zeros(p,m);
  sysh      = ss(a,b,c,d,-1);
  info.M    = 0;
  info.iter = 0;
  if nargout == 4
    xini = inistate(w,sysh);
  end
  return
end

% Initial approximation
if ~isempty(sys0)
  if ~isa(sys0,'ss')
    warning('OPT.SYS0 not an SS object. Ignored.')
    sys0 = [];
  else
    [pp,mm] = size(sys0); 
    nn = size(sys0,'order');
    if (mm ~= m) | (pp ~= p) | (nn ~= n)
      warning('OPT.SYS0 invalid. Ignored.')
      sys0 = [];
    end
  end
end
% If sys0 = [], replace it by a subspace estimate, otherwise TLS is used.
  
% Options for STLS
opt = [];
opt.epsrel  = epsrel;
opt.epsabs  = epsabs;
opt.epsgrad = epsgrad;
opt.maxiter = maxiter;
opt.disp    = disp;

% More constants
siso = (m == 1) & (p == 1);
l1   = l + 1; % number of block columns in the Hankel matrix

% Treat separately the general and the output only cases
if m > 0
  % Structure specification
  ne = length(exct);
  nn = nw - ne;
  noisy = 1:nw; noisy(exct) = []; % indices of noisy variables
  if ne == 0
    struct = [2 nw*l1 nw];
    h = blkhankel(w,l1)';    
  else % Define a block of exact variables
    struct = [4 ne*l1 ne; 2 nn*l1 nn]; 
    h = [blkhankel(w(:,exct,:),l1); blkhankel(w(:,noisy,:),l1)]';
  end
  if k > 1
    struct.a = struct;
    struct.k = k;
  end
  
  % Call the STLS solver
  x0 = sys2x(sys0, exct, noisy);
  [x,info] = stls(h(:,1:end-p), h(:,end-p+1:end), struct, x0, opt);
  info.M = sqrt(info.fmin);
  info = rmfield(info,'fmin');
  
  % Find the corrected trajectory
  if nargout > 2
    dp = corr(x, h, struct);
    clear h
    if k == 1
      dwn = reshape(dp,nn,T)';
    else % k > 1
      dwn = shiftdim(reshape(dp,nn,k,T),2);
    end
    w(:,noisy,:) = w(:,noisy,:) - dwn;
    clear dp dwn
  else % WH not needed
    clear h
  end

  % Find an I/S/O representation of the model
  R  = [x', -eye(p)]; % kernel representation of SYSH
  if ne > 0
    R3e = reshape(R(:,1:ne*l1),p,ne,l1);
    R3n = reshape(R(:,l1*ne+1:end),p,nn,l1);
    R3  = [R3e R3n];
    % Restore the original ordering
    R3(:,[exct noisy],:) = R3; 
    clear R3n R3e
  else % ne == 0
    R3 = reshape(R,p,nw,l1);
  end
  % Extract P and Q
  Q3 =  R3(:,1:m,:);   
  P3 = -R3(:,m+1:nw,:);
  a  = zeros(n); a(p+1:end,1:n-p) = eye(n-p);
  c  = [ zeros(p,n-p) eye(p) ]; 
  d  = Q3(:,:,l1);
  b  = zeros(n,m);
  ind_j = (n-p+1):n;
  for i = 1:l
    ind_i = (i-1)*p+1:i*p;
    P3i = P3(:,:,i);
    a(ind_i,ind_j) = - P3i;
    b(ind_i,:) = Q3(:,:,i) - P3i*d;
  end
  sysh = ss(a,b,c,d,-1);
  
  % Find the initial condition
  if nargout > 3
    xini = inistate(w,sysh);
  end
  
else % output only case (generically can not have exact variables)
  
  % Structure specification
  struct = [2 p*l1 p] ; 
  if k > 1
    struct.a = struct;
    struct.k = k;
  end
  
  % Call the stls solver
  h  = blkhankel(w,l1)';
  x0 = sys2x(sys0);
  [xh,info] = stls(h(:,1:end-p), h(:,end-p+1:end), struct, x0, opt);
  info.M = sqrt(info.fmin);
  info = rmfield(info,'fmin');
  
  % Find an input/state/output representation of the model
  O    = null([xh' -eye(p)]);
  ch   = O(1:p,:);
  ah   = O(1:end-p,:) \ O(p+1:end,:);
  sysh = ss(ah,[],ch,[],-1);      
  
  % Find the corrected trajectory
  if nargout > 2
    dp = corr(xh, h, struct);
    clear h
    uh = [];
    if k == 1
      dw = reshape(dp,p,T)';
    else % k > 1
      dw = shiftdim(reshape(dp,p,k,T),2);    
    end
    w = w - dw;
    clear dp dw
  end
  
  % Compute initial conditions
  if nargout > 3    
    xini = zeros(n,k);
    for i = 1:k
      tmp  = w(1:l1,:,i)'; 
      xini(:,i) = O \ tmp(:);
    end
  end
  
end

% ---------------------------------------------------------
% To do:
% 1. Include CORR in STLS. 
% 2. Native C-implementation of STLSIDENT and MISFIT.
% ---------------------------------------------------------