sgvm_getZdr

SGVM_GETZDR calculate driving point impedances
   ZDR = SGVM_GETZDR(MPC, FULL)

   Calculate the driving point impedances for the admittance matrix
   formed from the MATPOWER case, MPC. If FULL is true  the full inverse
   of the Ybus is calculated, otherwise a PCA version of the Ybus is used.

   Input
       MPC     matpower case
       FULL    1 for full inversion of Ybus,
               0 for PCA inversion
   Output
       ZDR     nb x 1 vector of complex driving point impedances

0001 function Zdr = sgvm_getZdr(mpc, full)
0002 %SGVM_GETZDR calculate driving point impedances
0003 %   ZDR = SGVM_GETZDR(MPC, FULL)
0004 %
0005 %   Calculate the driving point impedances for the admittance matrix
0006 %   formed from the MATPOWER case, MPC. If FULL is true  the full inverse
0007 %   of the Ybus is calculated, otherwise a PCA version of the Ybus is used.
0008 %
0009 %   Input
0010 %       MPC     matpower case
0011 %       FULL    1 for full inversion of Ybus,
0012 %               0 for PCA inversion
0013 %   Output
0014 %       ZDR     nb x 1 vector of complex driving point impedances
0015 
0016 %   SynGrid
0017 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0018 %   by Eran Schweitzer, Arizona State University
0019 %
0020 %   This file is part of SynGrid.
0021 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0022 
0023 if nargin < 2
0024     full = 0;
0025 end
0026 
0027 nb = size(mpc.bus,1);
0028 if ~all(mpc.bus(:,1) == (1:nb)')
0029     error(['sgvm_getZdr: bus numbers must be in internal order',...
0030  ' (consecutive starting at 1), consider using ext2int'])
0031 end
0032 
0033 Ybus = makeYbus(mpc);
0034 if full
0035     Zdr = diag(Ybus\speye(nb));
0036 else
0037     npca = min(100, round(0.1*nb));
0038     if have_fcn('octave')
0039         [U,S,V] = svds(Ybus, npca, 0);
0040     else
0041         [U,S,V] = svds(Ybus, npca, 'smallest');
0042     end
0043     if nb < 1e4
0044         Zdr = diag(V*(S\U'));
0045     else
0046         tmp = S\U';
0047         Zdr = zeros(nb,1);
0048         for k = 1:nb
0049             Zdr(k) = V(k,:)*tmp(:,k);
0050         end
0051     end
0052 end