stlsidentuy.m

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

function [sysh, info, uh, yh, xini] = stlsidentuy(u, y, l, oe, sys0, disp)
% STLSIDENTUY - Identification by structured total least squares.
%
% [sysh, info, uh, yh, xini] = stlsidentuy(u, y, l, oe, sys0, disp)
%
% U, Y   - input/output, for single time series TxM and TxP matrices.
%   For multiple time series, arrays with dimensions TxMxK and TxPxK.
%   If U = [], then an autonomous system is identified.
% L      - system lag.
% OE     - if nonzero OE identification (U is not approximated).
% SYS0   - initial approximation, an SS object of order L*P.
% DISP   - level of display [off | iter | notify | final].
% SYSH   - identified system, an SS object.
% 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.
% UH, YH - optimal approximating trajectory.
% XINI   - initial condition, under which (UH, YH) is obtained.

% Default options
opt.epsrel  = 1e-4;
opt.epsabs  = 1e-4;
opt.epsgrad = 1e-4;
opt.maxiter = 100;

if nargin == 6
  opt.disp = disp;
else % default
  opt.disp = 'notify';
end

% Constants
T    = size(y,1);
k    = size(y,3);
m    = size(u,2);
p    = size(y,2);
n    = l * p; 
l1   = l + 1;
siso = (m == 1) & (p == 1);

% Default OE = 0
if nargin < 4
  oe = 0;
elseif isempty(oe)
  oe = 0;
end

% Initial approximation
if nargin < 5
  sys0 = [];
end
x0 = sys2xuy(sys0);

if m > 0 % input-output identification
  
  % Structure specification
  if oe
    struct = [4 m*l1 m; 2 p*l1 p]; 
  else
    struct = [2 m*l1 m; 2 p*l1 p];
  end
  if k > 1
    struct.a = struct;
    struct.k = k;
  end  
  h = [blkhankel(u,l1)', blkhankel(y,l1)'];
    
  % Call the stls solver
  [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);
    if oe
      if k == 1
        dy = reshape(dp,p,T)';
      else
        dy = shiftdim(reshape(dp,p,k,T),2);
      end
      uh = u;
    else
      if k == 1
        du = reshape(dp(1:m*T),m,T)';
        dy = reshape(dp(m*T+1:end),p,T)';
      else
        du = shiftdim(reshape(dp(1:m*k*T),m,k,T),2);
        dy = shiftdim(reshape(dp(m*k*T+1:end),p,k,T),2);
      end
      uh = u - du;
      clear du
    end
    yh = y - dy;
    clear dp dy
  end
  clear h

  % Find the state space model
  Q  = x(1:l1*m,:)';
  P  = -x(l1*m+1:end,:)';
  a  = zeros(n);
  a(p+1:end,1:n-p) = eye(n-p);
  c  = [ zeros(p,n-p) eye(p) ]; 
  d  = Q(:,end-m+1:end); % = Q_{l+1}
  b  = zeros(n,m);
  ind_j = (n-p+1):n;
  for i = 1:l
    ind_i = (i-1)*p+1:i*p;
    Pi = P(:,ind_i);
    a(ind_i,ind_j) = - Pi;
    b(ind_i,:) = - Pi*d + Q(:,(i-1)*m+1:i*m);
  end
  sysh = ss(a,b,c,d,-1);

  % Find the initial condition
  if nargout > 4
    xini = inistate([uh yh], sysh);
  end
  
else % output only
  
  % Structure specification
  struct = [2 p*l1 p] ; 
  if k > 1
    struct.a = struct;
    struct.k = k;
  end  
  h = blkhankel(y,l1)';
  
  % Call the stls solver
  [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 state space model
  O   = null([x' -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(x, h, struct);
    clear h
    uh = [];
    if k == 1
      dy = reshape(dp,p,T)';
    else
      dy = shiftdim(reshape(dp,p,k,T),2);    
    end
    yh = y - dy;
    clear dp dy
  end
    
  % Find the initial condition
  if nargout > 4
    xini = zeros(n,k);
    for i = 1:k
      tmp  = y(1:l1,:,i)'; 
      xini(:,i) = O \ tmp(:);
    end
  end
  
end