misfituy.m

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

function [ M, uh, yh, xini ] = misfituy( u, y, sys, oe )
% MISFIT - Global Total Least Squares misfit.
% 
% [ M, uh, yh, xini ] = misfituy( u, y, sys, oe )
% 
% 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.
% SYS    - state-space system with M inputs, P outputs, 
%          and order N:=L*P, where L is an integer.
% OE     - if nonzero OE misfit (||Y - YH||_F).
% M      - misfit || [U Y] - [UH YH] ||_F.
% UH, YH - optimal approximating time series.
% XINI   - initial condition, under which (UH, YH) is obtained.

% Constants
T     = size(y,1); % # of samples
m     = size(u,2); % # of inputs
p     = size(y,2); % # of outputs
k     = size(y,3); % # of time series
n     = size(sys,'order');
l     = n/p;
l1    = l + 1;

if nargin < 4
  oe = 0;
end

% Structure specification
if m > 0
  if oe
    struct = [4 m*l1 m; 2 p*l1 p]; 
  else
    struct = [2 m*l1 m; 2 p*l1 p];
  end  
  h = [blkhankel(u,l1)', blkhankel(y,l1)'];
else % output only case
  struct = [2 p*l1 p] ; 
  h = blkhankel(y,l1)';
end
if k > 1
  struct.a = struct;
  struct.k = k;
end

% Convert (A,B,C,D) to a parameter X for STLS 
x = sys2xuy(sys);

% STLS cost function evaluation
if nargout > 1
  [M,dp] = stls_cost_fun(x,h,struct);
else
  M = stls_cost_fun(x,h,struct);
end
M = sqrt(M);
clear h

% Find the corrected trajectory
if nargout > 1
  if m > 0
    if oe
      if k == 1
        dy = reshape(dp,p,T)';
      else % k > 1
        dy = shiftdim(reshape(dp,p,k,T),2);    
      end
      uh = u;
    else % EIV
      if k == 1
        du = reshape(dp(1:m*T),m,T)';
        dy = reshape(dp(m*T+1:end),p,T)';
      else % k > 1
        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
    clear u
  else % output only case
    uh = [];
    if k == 1
      dy = reshape(dp,p,T)';
    else % k > 1
      dy = shiftdim(reshape(dp,p,k,T),2);    
    end
  end
  yh = y - dy;
  clear y dp dy
  
  % Find the initial condition
  if nargout > 2
    xini = inistate([uh yh],sys);
  else
    xini = NaN;
  end  
end