% Do A PBH test on system. Given a state space model SYS returns:
%
%  - Nr. of states
%  - Nr. of stable states
%  - Nr. of controllable states
%  - Nr. of unstabilizable states
%  - Nr. of observable states
%  - Nr. of undetectable states

function [nx, nsta, nctrb, nustab, nobsv, nudetb] = pbhtest(sys)

nx = length(sys.A);
eigvals = eig(sys.A);

% Count number of stable states
nsta = 0;
for i = 1:nx
  if real(eigvals(i)) < 0
    nsta = nsta + 1;
  end
end

% Check system controllability / stabilizability
Wc = ctrb(sys);
nctrb = rank(Wc);
nustab = 0;
if nctrb < nx
  % Is the system at least stabilizable?
  for i = 1:nx
    if real(eigvals(i)) >= 0
      % PBH test
      W = [(sys.A - eigvals(i) * eye(nx)), sys.B];
      if rank(W) < nx
        nustab = nustab + 1;
      end
    end
  end
end

% Check system observability / detectability
Wo = obsv(sys);
nobsv = rank(Wo);
nudetb = 0;
if nobsv < nx
  % is the system at least detectable?
  for i = 1:nx
    if real(eigvals(i)) >= 0
      % PBH test
      W = [(sys.A' - eigvals(i) * eye(nx)), sys.C'];
      if rank(W) < nx
        nudetb = nudetb + 1;
      end
    end
  end
end

end