MAKEPTDF Builds the DC PTDF matrix for a given choice of slack. H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK) 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. Examples: H = makePTDF(baseMVA, bus, branch); H = makePTDF(baseMVA, bus, branch, 1); slack = rand(size(bus, 1), 1); H = makePTDF(baseMVA, bus, branch, slack); See also MAKELODF.
0001 function H = makePTDF(baseMVA, bus, branch, slack) 0002 %MAKEPTDF Builds the DC PTDF matrix for a given choice of slack. 0003 % H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK) returns the DC PTDF 0004 % matrix for a given choice of slack. The matrix is nbr x nb, where 0005 % nbr is the number of branches and nb is the number of buses. The SLACK 0006 % can be a scalar (single slack bus) or an nb x 1 column vector of 0007 % weights specifying the proportion of the slack taken up at each bus. 0008 % If the SLACK is not specified the reference bus is used by default. 0009 % 0010 % Examples: 0011 % H = makePTDF(baseMVA, bus, branch); 0012 % H = makePTDF(baseMVA, bus, branch, 1); 0013 % slack = rand(size(bus, 1), 1); 0014 % H = makePTDF(baseMVA, bus, branch, slack); 0015 % 0016 % See also MAKELODF. 0017 0018 % For convenience, SLACK can also be an nb x nb matrix, where each 0019 % column specifies how the slack should be handled for injections 0020 % at that bus. 0021 0022 % MATPOWER 0023 % $Id: makePTDF.m 1635 2010-04-26 19:45:26Z ray $ 0024 % by Ray Zimmerman, PSERC Cornell 0025 % Copyright (c) 2006-2010 by Power System Engineering Research Center (PSERC) 0026 % 0027 % This file is part of MATPOWER. 0028 % See http://www.pserc.cornell.edu/matpower/ for more info. 0029 % 0030 % MATPOWER is free software: you can redistribute it and/or modify 0031 % it under the terms of the GNU General Public License as published 0032 % by the Free Software Foundation, either version 3 of the License, 0033 % or (at your option) any later version. 0034 % 0035 % MATPOWER is distributed in the hope that it will be useful, 0036 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0037 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0038 % GNU General Public License for more details. 0039 % 0040 % You should have received a copy of the GNU General Public License 0041 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0042 % 0043 % Additional permission under GNU GPL version 3 section 7 0044 % 0045 % If you modify MATPOWER, or any covered work, to interface with 0046 % other modules (such as MATLAB code and MEX-files) available in a 0047 % MATLAB(R) or comparable environment containing parts covered 0048 % under other licensing terms, the licensors of MATPOWER grant 0049 % you additional permission to convey the resulting work. 0050 0051 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0052 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0053 0054 %% use reference bus for slack by default 0055 if nargin < 4 0056 slack = find(bus(:, BUS_TYPE) == REF); 0057 slack = slack(1); 0058 end 0059 0060 %% set the slack bus to be used to compute initial PTDF 0061 if length(slack) == 1 0062 slack_bus = slack; 0063 else 0064 slack_bus = 1; %% use bus 1 for temp slack bus 0065 end 0066 0067 nb = size(bus, 1); 0068 nbr = size(branch, 1); 0069 noref = (2:nb)'; %% use bus 1 for voltage angle reference 0070 noslack = find((1:nb)' ~= slack_bus); 0071 0072 %% check that bus numbers are equal to indices to bus (one set of bus numbers) 0073 if any(bus(:, BUS_I) ~= (1:nb)') 0074 error('makePTDF: buses must be numbered consecutively in bus matrix') 0075 end 0076 0077 %% compute PTDF for single slack_bus 0078 [Bbus, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch); 0079 H = zeros(nbr, nb); 0080 H(:, noslack) = full(Bf(:, noref) / Bbus(noslack, noref)); 0081 %% = full(Bf(:, noref) * inv(Bbus(noslack, noref))); 0082 0083 %% distribute slack, if requested 0084 if length(slack) ~= 1 0085 if size(slack, 2) == 1 %% slack is a vector of weights 0086 slack = slack/sum(slack); %% normalize weights 0087 0088 %% conceptually, we want to do ... 0089 %% H = H * (eye(nb,nb) - slack * ones(1, nb)); 0090 %% ... we just do it more efficiently 0091 v = H * slack; 0092 for k = 1:nb 0093 H(:, k) = H(:, k) - v; 0094 end 0095 else 0096 H = H * slack; 0097 end 0098 end