Home > matpower7.1 > extras > syngrid > lib > sgvm_acptdf.m

sgvm_acptdf

PURPOSE ^

SGVM_ACPTDF Calculate line sensitivitis at a given operating point.

SYNOPSIS ^

function [Hf, Ht, threshold] = sgvm_acptdf(mpc, bidx, pca)

DESCRIPTION ^

SGVM_ACPTDF Calculate line sensitivitis at a given operating point.
   [HF, HT, THRESHOLD] = SGVM_ACPTDF(MPC, BIDX, PCA)

   We want:
          [dPf; dQf]  = Hf*[dPbus; dQbus]  (1)
   Hf is a 2nl x 2(nb-1) sensitivity matrix:
          | dPf/dPbus     dPf/dQbus |
    Hf =  |                         |
          | dQf/dPbus     dQf/dQbus |
   now,
          [dPbus; dQbus] = J*[dVa ; dVm]
   and,
          [dPf; dQf] = Kf*[dVa ; dVm]
   where,
          | dPbus/dVa     dPbus/dVm |
     J =  |                         |
          | dQbus/dVa     dQbus/dVm |
   and,
          | dPf/dVa       dPf/dVm |
    Kf =  |                       |
          | dQf/dVa       dQf/dVm |
   Then (1) can be expressed as
     Kf = Hf*J -----> Hf = Kf*J^-1
   Note that as presented here, these calculations are slack bus
   dependent.

   Inputs
     MPC  - MATPOWER case (should be in internal indexing mode)
     BIDX - branch ids to use (default is all branches)
     PCA  - Number of principle components to use when inverting the Jacobian matrix.
              0 - use the full jacobian
              pca > 0 - use the 'pca' smallest singular values of J to estimate J^-1

   Outputs
     HF        - line sensitivities matrix using from end powers.
     HT        - line sensitivities matrix using to end powers.
     THRESHOLD - values with absolute value below THRESHOLD are set to 0.
                 The THRESHOLD is 5% of the maximum magnitude in HF or HF and HT
                 if both are returned.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Hf, Ht, threshold] = sgvm_acptdf(mpc, bidx, pca)
0002 %SGVM_ACPTDF Calculate line sensitivitis at a given operating point.
0003 %   [HF, HT, THRESHOLD] = SGVM_ACPTDF(MPC, BIDX, PCA)
0004 %
0005 %   We want:
0006 %          [dPf; dQf]  = Hf*[dPbus; dQbus]  (1)
0007 %   Hf is a 2nl x 2(nb-1) sensitivity matrix:
0008 %          | dPf/dPbus     dPf/dQbus |
0009 %    Hf =  |                         |
0010 %          | dQf/dPbus     dQf/dQbus |
0011 %   now,
0012 %          [dPbus; dQbus] = J*[dVa ; dVm]
0013 %   and,
0014 %          [dPf; dQf] = Kf*[dVa ; dVm]
0015 %   where,
0016 %          | dPbus/dVa     dPbus/dVm |
0017 %     J =  |                         |
0018 %          | dQbus/dVa     dQbus/dVm |
0019 %   and,
0020 %          | dPf/dVa       dPf/dVm |
0021 %    Kf =  |                       |
0022 %          | dQf/dVa       dQf/dVm |
0023 %   Then (1) can be expressed as
0024 %     Kf = Hf*J -----> Hf = Kf*J^-1
0025 %   Note that as presented here, these calculations are slack bus
0026 %   dependent.
0027 %
0028 %   Inputs
0029 %     MPC  - MATPOWER case (should be in internal indexing mode)
0030 %     BIDX - branch ids to use (default is all branches)
0031 %     PCA  - Number of principle components to use when inverting the Jacobian matrix.
0032 %              0 - use the full jacobian
0033 %              pca > 0 - use the 'pca' smallest singular values of J to estimate J^-1
0034 %
0035 %   Outputs
0036 %     HF        - line sensitivities matrix using from end powers.
0037 %     HT        - line sensitivities matrix using to end powers.
0038 %     THRESHOLD - values with absolute value below THRESHOLD are set to 0.
0039 %                 The THRESHOLD is 5% of the maximum magnitude in HF or HF and HT
0040 %                 if both are returned.
0041 
0042 %   SynGrid
0043 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0044 %   by Eran Schweitzer, Arizona State University
0045 %
0046 %   This file is part of SynGrid.
0047 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0048 
0049 %% define indices
0050 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0051     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0052 nb = size(mpc.bus,1);
0053 nl = size(mpc.branch,1);
0054 
0055 if any(mpc.bus(:,BUS_I) ~= (1:nb)')
0056     error('sgvm_acptdf: bus matrix needs to have consecutive bus numbering.')
0057 end
0058 
0059 if nargin < 3
0060     pca = 0;
0061     if nargin < 2
0062         bidx = (1:nl).';
0063     else
0064         nl = length(bidx);
0065     end
0066 end
0067 %% calculate derivatives at operating point
0068 [Ybus, Yf, Yt] = makeYbus(mpc);
0069 V = mpc.bus(:,VM).*exp(1i*mpc.bus(:,VA)*pi/180);
0070 [dsf_dva, dsf_dvm, dst_dva, dst_dvm] = dSbr_dV(mpc.branch, Yf, Yt, V);
0071 [dsbus_dva, dsbus_dvm] = dSbus_dV(Ybus,V);
0072 
0073 %% form matrices and solve
0074 nonslack = find(mpc.bus(:,BUS_TYPE) ~= REF);
0075 % make Jacobian
0076 J = [real(dsbus_dva(nonslack,nonslack)), real(dsbus_dvm(nonslack,nonslack));
0077      imag(dsbus_dva(nonslack,nonslack)), imag(dsbus_dvm(nonslack,nonslack))];
0078 % make branch Jacobian
0079 Kf = [real(dsf_dva(bidx,nonslack)), real(dsf_dvm(bidx,nonslack));
0080       imag(dsf_dva(bidx,nonslack)), imag(dsf_dvm(bidx,nonslack))];
0081 Kt = [real(dst_dva(bidx,nonslack)), real(dst_dvm(bidx,nonslack));
0082       imag(dst_dva(bidx,nonslack)), imag(dst_dvm(bidx,nonslack))];
0083 
0084 % solve for sensitivities.
0085 Hf = sparse(2*nl, 2*nb);
0086 if pca > 0
0087     [U,S,V] = svds(J,pca,'smallest');
0088     Hf(:,[nonslack; nb+nonslack]) = Kf * (V*(S\U'));
0089 else
0090     Hf(:,[nonslack; nb+nonslack]) = Kf / J;
0091 end
0092 if nargout > 1
0093     Ht = sparse(2*nl, 2*nb);
0094     if pca
0095         Ht(:,[nonslack; nb+nonslack]) = Kt * (V*(S\U'));
0096     else
0097         Ht(:,[nonslack; nb+nonslack]) = Kt / J;
0098     end
0099 end
0100 
0101 % threshold
0102 if nargout > 1
0103 %   threshold = full(0.05*max(abs([Hf(:);Ht(:)])));
0104 %     threshold = min(0.01, threshold);
0105     threshold = 0.05*max(max(abs(Hf), abs(Ht)),[], 2);
0106     if have_fcn('octave')
0107         Ht(full(abs(Ht)) < full(threshold)) = 0;
0108     else
0109         Ht(abs(Ht) < threshold) = 0;
0110     end
0111 else
0112 %   threshold = full(0.05*max(abs(Hf(:))));
0113 %   threshold = min(0.01, threshold);
0114     threshold = 0.05*max(abs(Hf),[], 2);
0115 end
0116 if have_fcn('octave')
0117     Hf(full(abs(Hf)) < full(threshold)) = 0;
0118 else
0119     Hf(abs(Hf) < threshold) = 0;
0120 end

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005