makePTDF

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,v 1.9 2010/04/26 19:45:25 ray Exp $
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