MAKEPTDF Builds the DC PTDF matrix for a given choice of slack.
H = MAKEPTDF(MPC)
H = MAKEPTDF(MPC, SLACK)
H = MAKEPTDF(MPC, SLACK, TXFR)
H = MAKEPTDF(MPC, SLACK, BUS_IDX)
H = MAKEPTDF(BASEMVA, BUS, BRANCH)
H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK)
H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK, TXFR)
H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK, BUS_IDX)
Returns the DC PTDF matrix for a given choice of slack. The matrix is
nbr x nb, where nbr is the number of branches and nb is the number of
buses. The SLACK can be a scalar (single slack bus) or an nb x 1 column
vector of weights specifying the proportion of the slack taken up at each
bus. If the SLACK is not specified the reference bus is used by default.
Bus numbers must be consecutive beginning at 1 (i.e. internal ordering).
If TXFR is supplied it must be a matrix (or vector) with nb rows whose
columns each sum to zero, where each column defines a specific (slack
independent) transfer. E.g. if k-th transfer is from bus i to bus j,
TXFR(i, k) = 1 and TXFR(j, k) = -1. In this case H has the same number
of columns as TXFR.
If BUS_IDX is supplied, it contains a column vector of bus indices.
The columns of H correspond to these indices, but they are computed
individually rather than computing the full PTDF matrix and selecting
the desired columns.
Examples:
H = makePTDF(mpc);
H = makePTDF(baseMVA, bus, branch, 1);
slack = rand(size(bus, 1), 1);
H = makePTDF(mpc, slack);
% for transfer from bus i to bus j
txfr = zeros(nb, 1); txfr(i) = 1; txfr(j) = -1;
H = makePTDF(mpc, slack, txfr);
% for buses i and j only
H = makePTDF(mpc, slack, [i;j]);
See also MAKELODF.0001 function H = makePTDF(baseMVA, bus, branch, slack, bus_idx) 0002 %MAKEPTDF Builds the DC PTDF matrix for a given choice of slack. 0003 % H = MAKEPTDF(MPC) 0004 % H = MAKEPTDF(MPC, SLACK) 0005 % H = MAKEPTDF(MPC, SLACK, TXFR) 0006 % H = MAKEPTDF(MPC, SLACK, BUS_IDX) 0007 % H = MAKEPTDF(BASEMVA, BUS, BRANCH) 0008 % H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK) 0009 % H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK, TXFR) 0010 % H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK, BUS_IDX) 0011 % 0012 % Returns the DC PTDF matrix for a given choice of slack. The matrix is 0013 % nbr x nb, where nbr is the number of branches and nb is the number of 0014 % buses. The SLACK can be a scalar (single slack bus) or an nb x 1 column 0015 % vector of weights specifying the proportion of the slack taken up at each 0016 % bus. If the SLACK is not specified the reference bus is used by default. 0017 % Bus numbers must be consecutive beginning at 1 (i.e. internal ordering). 0018 % 0019 % If TXFR is supplied it must be a matrix (or vector) with nb rows whose 0020 % columns each sum to zero, where each column defines a specific (slack 0021 % independent) transfer. E.g. if k-th transfer is from bus i to bus j, 0022 % TXFR(i, k) = 1 and TXFR(j, k) = -1. In this case H has the same number 0023 % of columns as TXFR. 0024 % 0025 % If BUS_IDX is supplied, it contains a column vector of bus indices. 0026 % The columns of H correspond to these indices, but they are computed 0027 % individually rather than computing the full PTDF matrix and selecting 0028 % the desired columns. 0029 % 0030 % Examples: 0031 % H = makePTDF(mpc); 0032 % H = makePTDF(baseMVA, bus, branch, 1); 0033 % slack = rand(size(bus, 1), 1); 0034 % H = makePTDF(mpc, slack); 0035 % 0036 % % for transfer from bus i to bus j 0037 % txfr = zeros(nb, 1); txfr(i) = 1; txfr(j) = -1; 0038 % H = makePTDF(mpc, slack, txfr); 0039 % 0040 % % for buses i and j only 0041 % H = makePTDF(mpc, slack, [i;j]); 0042 % 0043 % See also MAKELODF. 0044 0045 % For convenience, SLACK can also be an nb x nb matrix, where each 0046 % column specifies how the slack should be handled for injections 0047 % at that bus. This option only applies when computing the full 0048 % PTDF matrix (i.e. when TXFR and BUS_IDX are not provided.) 0049 0050 % MATPOWER 0051 % Copyright (c) 2006-2020, Power Systems Engineering Research Center (PSERC) 0052 % by Ray Zimmerman, PSERC Cornell 0053 % 0054 % This file is part of MATPOWER. 0055 % Covered by the 3-clause BSD License (see LICENSE file for details). 0056 % See https://matpower.org for more info. 0057 0058 %% extract from MPC if necessary 0059 if isstruct(baseMVA) 0060 mpc = baseMVA; 0061 if nargin < 3 0062 bus_idx = []; 0063 if nargin < 2 0064 slack = []; 0065 else 0066 slack = bus; 0067 end 0068 else 0069 slack = bus; 0070 bus_idx = branch; 0071 end 0072 baseMVA = mpc.baseMVA; 0073 bus = mpc.bus; 0074 branch = mpc.branch; 0075 else 0076 if nargin < 5 0077 bus_idx = []; 0078 if nargin < 4 0079 slack = []; 0080 end 0081 end 0082 end 0083 0084 %% define named indices into bus matrix 0085 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0086 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0087 0088 %% use reference bus for slack by default 0089 if isempty(slack) 0090 slack = find(bus(:, BUS_TYPE) == REF); 0091 slack = slack(1); 0092 end 0093 0094 %% compute full PTDF? 0095 nb = size(bus, 1); 0096 nbr = size(branch, 1); 0097 txfr = 0; %% default assumes not just for specific transfers 0098 if ~isempty(bus_idx) 0099 compute_full_H = 0; 0100 if size(bus_idx, 1) == nb && sum(sum(bus_idx)) == 0 0101 txfr = 1; %% cols of H for specific (slack independent) transfers 0102 dP = bus_idx; 0103 end 0104 else 0105 compute_full_H = 1; 0106 end 0107 0108 %% set the slack bus to be used to compute initial PTDF 0109 if length(slack) == 1 0110 slack_bus = slack; 0111 else 0112 slack_bus = 1; %% use bus 1 for temp slack bus 0113 end 0114 0115 noref = (2:nb)'; %% use bus 1 for voltage angle reference 0116 noslack = find((1:nb)' ~= slack_bus); 0117 0118 %% check that bus numbers are equal to indices to bus (one set of bus numbers) 0119 if any(bus(:, BUS_I) ~= (1:nb)') 0120 error('makePTDF: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering') 0121 end 0122 0123 %%----- compute PTDF for single slack_bus ----- 0124 [Bbus, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch); 0125 0126 %% set up injections/transfers 0127 if compute_full_H %% full H for all columns 0128 % %% old (slower) method 0129 % H = zeros(nbr, nb); 0130 % H(:, noslack) = full(Bf(:, noref) / Bbus(noslack, noref)); 0131 % %% = full(Bf(:, noref) * inv(Bbus(noslack, noref))); 0132 nbi = nb; %% number of buses of interest (i.e. all of them) 0133 dP = speye(nb, nb); 0134 else %% compute only for specific columns 0135 if txfr %% for specific transfers 0136 nbi = size(dP, 2); %% number of transfers 0137 else %% for subset of columns of H for specific buses 0138 %% if necessary add missing slacks to bus index list 0139 nbi0 = length(bus_idx); %% number of original buses of interest 0140 bidx = bus_idx; 0141 if length(slack) ~= 1 && size(slack, 2) == 1 0142 slacks = find(slack); %% get all slack buses 0143 k = find(~ismember(slacks, bus_idx)); %% find slacks not in bus_idx 0144 if ~isempty(k) 0145 bidx = [bus_idx; slacks(k)]; 0146 end 0147 end 0148 nbi = size(bidx, 1); %% number of buses of interest 0149 0150 %% define the equivalent transfer, each column is a single transfer 0151 dP = accumarray([bidx (1:nbi)'], 1, [nb, nbi]); 0152 end 0153 end 0154 0155 %% solve for change in voltage angles 0156 dTheta = zeros(nb, nbi); 0157 dTheta(noref, :) = Bbus(noslack, noref) \ dP(noslack, :); 0158 0159 %% compute corresponding change in branch flows 0160 H = Bf * dTheta; 0161 0162 %%----- distribute slack, if requested ----- 0163 if length(slack) ~= 1 && ~txfr 0164 if size(slack, 2) == 1 %% slack is a vector of weights 0165 slack = slack/sum(slack); %% normalize weights 0166 0167 %% conceptually, we want to do ... 0168 %% H = H * (eye(nb,nb) - slack * ones(1, nb)); 0169 %% ... we just do it more efficiently 0170 if compute_full_H 0171 v = H * slack; 0172 for k = 1:nb 0173 H(:, k) = H(:, k) - v; 0174 end 0175 else 0176 v = H * slack(bidx); 0177 for k = 1:nbi 0178 H(:, k) = H(:, k) - v; 0179 end 0180 H = H(:, 1:nbi0); %% remove temp cols added for missing slacks 0181 end 0182 else 0183 H = H * slack; 0184 end 0185 end