isobservable

ISOBSERVABLE  Test for observability.
   returns 1 if the system is observable, 0 otherwise.
   created by Rui Bo on Jan 9, 2010

0001 function TorF = isobservable(H, pv, pq)
0002 %ISOBSERVABLE  Test for observability.
0003 %   returns 1 if the system is observable, 0 otherwise.
0004 %   created by Rui Bo on Jan 9, 2010
0005 
0006 %   MATPOWER
0007 %   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
0008 %   by Rui Bo
0009 %
0010 %   This file is part of MATPOWER/mx-se.
0011 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0012 %   See https://github.com/MATPOWER/mx-se/ for more info.
0013 
0014 %% options
0015 tol     = 1e-5; % mpopt.pf.tol;
0016 check_reason = 1;    % check reason for system being not observable
0017                     % 0: no check
0018                     % 1: check (NOTE: may be time consuming due to svd calculation)
0019 
0020 %% test if H is full rank
0021 [m, n]  = size(H);
0022 r       = rank(H);
0023 if r < min(m, n)
0024     TorF = 0;
0025 else
0026     TorF = 1;
0027 end
0028 
0029 %% look for reasons for system being not observable
0030 if check_reason && ~TorF
0031     %% look for variables not being observed
0032     idx_trivialColumns = [];
0033     varNames = {};
0034     for j = 1:n
0035         normJ = norm(H(:, j), inf);
0036         if normJ < tol % found a zero column
0037             idx_trivialColumns = [idx_trivialColumns j];
0038             varName = getVarName(j, pv, pq);
0039             varNames{length(idx_trivialColumns)} = varName;
0040         end
0041     end
0042 
0043     if ~isempty(idx_trivialColumns) % found zero-valued column vector
0044         fprintf('Warning: The following variables are not observable since they are not related with any measurement!');
0045         varNames
0046         idx_trivialColumns
0047     else % no zero-valued column vector
0048         %% look for dependent column vectors
0049         for j = 1:n
0050             rr = rank(H(:, 1:j));
0051             if rr ~= j % found dependent column vector
0052                 %% look for linearly depedent vector
0053                 colJ = H(:, j); % j(th) column of H
0054                 varJName = getVarName(j, pv, pq);
0055                 for k = 1:j-1
0056                     colK = H(:, k);
0057                     if rank([colK colJ]) < 2 % k(th) column vector is linearly dependent of j(th) column vector
0058                         varKName = getVarName(k, pv, pq);
0059                         fprintf('Warning: %d(th) column vector (w.r.t. %s) of H is linearly dependent of %d(th) column vector (w.r.t. %s)!\n', j, varJName, k, varKName);
0060                         return;
0061                     end
0062                 end
0063             end
0064         end
0065     fprintf('Warning: No specific reason was found for system being not observable.\n');
0066     end
0067 end