sys2x.m

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

function x = sys2x(sys, exct, noisy)
% SYS2X - 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)];
    % Reorder as R = [Q0 P0 ... Ql Pl]
    R([1:2:end,2:2:end]) = R;
    % A tensor form
    R3 = reshape(R,1,2,l1);
  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;
    R3 = [reshape(Q,p,m,l1) -reshape(P,p,p,l1)];
  end
  % Reorder first exact then noisy
  ne = length(exct); nn = length(noisy);
  R  = [reshape(R3(:,exct,:),p,l1*ne) reshape(R3(:,noisy,:),p,l1*nn)]';
else % Output only (generically no exact variables)
  % 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,:);