makePTDF

MAKEPTDF   Builds the DC PTDF matrix for a given choice of slack.
   H = MAKEPTDF(MPC)
   H = MAKEPTDF(MPC, SLACK)
   H = MAKEPTDF(BASEMVA, BUS, BRANCH)
   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.
   Bus numbers must be consecutive beginning at 1 (i.e. internal ordering).

   Examples:
       H = makePTDF(mpc);
       H = makePTDF(baseMVA, bus, branch, 1);
       slack = rand(size(bus, 1), 1);
       H = makePTDF(mpc, 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(MPC)
0004 %   H = MAKEPTDF(MPC, SLACK)
0005 %   H = MAKEPTDF(BASEMVA, BUS, BRANCH)
0006 %   H = MAKEPTDF(BASEMVA, BUS, BRANCH, SLACK)
0007 %
0008 %   Returns the DC PTDF matrix for a given choice of slack. The matrix is
0009 %   nbr x nb, where nbr is the number of branches and nb is the number of
0010 %   buses. The SLACK can be a scalar (single slack bus) or an nb x 1 column
0011 %   vector of weights specifying the proportion of the slack taken up at each
0012 %   bus. If the SLACK is not specified the reference bus is used by default.
0013 %   Bus numbers must be consecutive beginning at 1 (i.e. internal ordering).
0014 %
0015 %   Examples:
0016 %       H = makePTDF(mpc);
0017 %       H = makePTDF(baseMVA, bus, branch, 1);
0018 %       slack = rand(size(bus, 1), 1);
0019 %       H = makePTDF(mpc, slack);
0020 %
0021 %   See also MAKELODF.
0022 
0023 %   For convenience, SLACK can also be an nb x nb matrix, where each
0024 %   column specifies how the slack should be handled for injections
0025 %   at that bus.
0026 
0027 %   MATPOWER
0028 %   Copyright (c) 2006-2016, Power Systems Engineering Research Center (PSERC)
0029 %   by Ray Zimmerman, PSERC Cornell
0030 %
0031 %   This file is part of MATPOWER.
0032 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0033 %   See https://matpower.org for more info.
0034 
0035 %% extract from MPC if necessary
0036 if nargin < 3
0037     mpc     = baseMVA;
0038     if nargin == 2
0039         slack = bus;
0040     end
0041     baseMVA = mpc.baseMVA;
0042     bus     = mpc.bus;
0043     branch  = mpc.branch;
0044 end
0045 
0046 %% define named indices into bus matrix
0047 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0048     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0049 
0050 %% use reference bus for slack by default
0051 if nargin == 1 || nargin == 3
0052     slack = find(bus(:, BUS_TYPE) == REF);
0053     slack = slack(1);
0054 end
0055 
0056 %% set the slack bus to be used to compute initial PTDF
0057 if length(slack) == 1
0058     slack_bus = slack;
0059 else
0060     slack_bus = 1;      %% use bus 1 for temp slack bus
0061 end
0062 
0063 nb = size(bus, 1);
0064 nbr = size(branch, 1);
0065 noref   = (2:nb)';      %% use bus 1 for voltage angle reference
0066 noslack = find((1:nb)' ~= slack_bus);
0067 
0068 %% check that bus numbers are equal to indices to bus (one set of bus numbers)
0069 if any(bus(:, BUS_I) ~= (1:nb)')
0070     error('makePTDF: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering')
0071 end
0072 
0073 %% compute PTDF for single slack_bus
0074 [Bbus, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0075 H = zeros(nbr, nb);
0076 H(:, noslack) = full(Bf(:, noref) / Bbus(noslack, noref));
0077         %%    = full(Bf(:, noref) * inv(Bbus(noslack, noref)));
0078 
0079 %% distribute slack, if requested
0080 if length(slack) ~= 1
0081     if size(slack, 2) == 1  %% slack is a vector of weights
0082         slack = slack/sum(slack);   %% normalize weights
0083         
0084         %% conceptually, we want to do ...
0085         %%    H = H * (eye(nb,nb) - slack * ones(1, nb));
0086         %% ... we just do it more efficiently
0087         v = H * slack;
0088         for k = 1:nb
0089             H(:, k) = H(:, k) - v;
0090         end
0091     else
0092         H = H * slack;
0093     end
0094 end