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.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