sys2xuy.m

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

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

if ~isa(sys,'ss')
  x = []; % Use default for the stls solver
  return
end

[p,m] = size(sys);
n     = size(sys,'order');
l     = n / p;
l1    = l + 1;

if (m > 0)
  if (p == 1) & (m == 1)
    [Q,P] = tfdata(sys,'v');
    R     = [fliplr(Q), -fliplr(P)]';
  else % MIMO
    % Compute the observability matrix
    O = zeros(l1*p,n);
    O(1:p,:) = sys.c;
    for i = 2:l1
      O((i-1)*p+1:i*p,:) = O((i-2)*p+1:(i-1)*p,:) * sys.a;
    end
    P = null(O')';
    % Form the lower bock-triangular Toeplitz matrix T whose first 
    % block column is F = [ D; CB; CAB; ...; CA^{l-1}B ]
    F = [ sys.d; O(1:end-p,:) * sys.b ];
    T = zeros(l1*p,l1*m);
    for i = 1:l1
      T((i-1)*p+1:end,(i-1)*m+1:i*m) = F(1:(l1+1-i)*p,:);
    end
    Q = P * T;
    R = [Q -P]';
  end
else % Output only
  % Compute the observability matrix
  O = zeros(l1*p,n);
  O(1:p,:) = sys.c;
  for t = 2:l1
    O((t-1)*p+1:t*p,:) = O((t-2)*p+1:(t-1)*p,:)*sys.a;
  end
  R = null(O');
end

% Normalize with last block -I, in order to get X
x = -R(1:end-p,:)/R(end-p+1:end,:);