misfit.m

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

function [ M, w, xini ] = misfit( w, sys, exct )
% MISFIT - Global Total Least Squares misfit ||W - Wh||_F. 
% 
% [ M, wh, xini ] = misfit( w, sys, exct )
% 
% Inputs:
% W    - given time series (#samples x #variables x #time series).
% SYS  - state-space system with M inputs, P:=size(W,2)-M outputs, 
%        and order N:=L*P, where L is an integer.
% EXCT - vector of indices for exact variables (default []).
%
% Outputs:
% M    - misfit ||W - WH||_F.
% WH   - optimal approximating time series.
% XINI - initial condition, under which WH is obtained.

% 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

% Constants
T     = size(w,1); % # of samples
nw    = size(w,2); % # of variables
k     = size(w,3); % # of time series
[p,m] = size(sys);
n     = size(sys,'order');
l     = n/p;
l1    = l + 1;

% Check for parameters out of range
if nargin < 2
  error('Not enough input arguments. (W and SYS should be given.)')
end
if (~isreal(w))
  error('Complex data not supported.')
end
if nargin == 2 | isempty(exct)
  exct = [];
elseif any(1 > exct | nw < exct)
  error('An index for exact variable is out of range.')
elseif m == 0
  error('The output only case can not have exact variables.')
end

% Check SYS
if ~isa(sys,'ss')
  error('SYS not an SS object.')
end
if (m+p ~= nw)
  error('#inputs + #output of SYS not equal to size(W,2).')
end
if (ceil(l)-l ~=  0)
  error('SYS not of order that is a multiple of #outputs.')
end

% Structure specification
noisy = 1:nw; % indices of noisy variables
if m > 0
  ne = length(exct);
  nn = nw - ne;
  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
else % Output only case (generically can not have exact variables)
  struct = [2 p*l1 p];
  h = blkhankel(w,l1)';     
end
if k > 1
  struct.a = struct;
  struct.k = k;
end  

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

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

% Find the corrected trajectory
if nargout > 1
  if m > 0 % Input-Output case
    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;
  else % Output only case
    if k == 1
      dw = reshape(dp,p,T)';
    else % k > 1
      dw = shiftdim(reshape(dp,p,k,T),2);    
    end
    w = w - dw;
  end
  % Find the initial condition
  if nargout > 2
    xini = inistate(w,sys);
  end  
end