% Generates problem for system identification
%
% x_{t+1}  = A*x_t + B*u_t
% bx_t     = C*x_t + D*u_t
%
% where A is r by r, B is r by p, C is m by r, D is m by p
% the input u is p by N+1 and the output bx is m by N+1
%

  N = 1299; r = 10; m = 20;
  p = 5;
  
  A = randn(r,r);
  A = A/norm(A); % normalize A
  B = randn(r,p);
  B = B/norm(B); % normalize B
  C = randn(m,r);
  C = C/norm(C);
  D = randn(m,p);
  D = D/norm(D);
  u = randn(p,N+1);
  
  x = randn(r,1); % randomly generate state variable
  
  bx = zeros(m,N+1); % collection of all outputs
  
  for i = 1:N+1;
    xnew = A*x + B*u(:,i);
    bx(:,i) = C*x + D*u(:,i);
    x = xnew;
  end
  
  bx = bx + 5e-2*randn(m,N+1); % add noise
   % the r above are real system order. the r below is the number of blocks
   % required
  r = 2*r + 1;
  U = Hrx(u,r);
  U = null(U);
  q = size(U,2);
  
%   save sysP8 m N p r U bx