makeJac

MAKEJAC  Forms the power flow Jacobian.
   J = MAKEJAC(MPC)
   J = MAKEJAC(MPC, FULLJAC)
   J = MAKEJAC(BASEMVA, BUS, BRANCH, GEN)
   J = MAKEJAC(BASEMVA, BUS, BRANCH, GEN, FULLJAC)
   [J, YBUS, YF, YT] = MAKEJAC(MPC)

   Returns the power flow Jacobian and, optionally, the system admittance
   matrices. Inputs can be a MATPOWER case struct or individual BASEMVA,
   BUS, BRANCH and GEN values. Bus numbers must be consecutive beginning
   at 1 (i.e. internal ordering). If the FULLJAC argument is present and
   true, it returns the full Jacobian (sensitivities of all bus injections
   w.r.t all voltage angles/magnitudes) as opposed to the reduced version
   used in the Newton power flow updates. The units for all quantities are
   in per unit with radians for voltage angles.

   Note: This function builds the Jacobian from scratch, rebuilding the
         YBUS matrix in the process. You probably don't want to use this
         in performance critical code.

   See also MAKEYBUS, EXT2INT

0001 function [J, Ybus, Yf, Yt] = makeJac(baseMVA, bus, branch, gen, fullJac)
0002 %MAKEJAC  Forms the power flow Jacobian.
0003 %   J = MAKEJAC(MPC)
0004 %   J = MAKEJAC(MPC, FULLJAC)
0005 %   J = MAKEJAC(BASEMVA, BUS, BRANCH, GEN)
0006 %   J = MAKEJAC(BASEMVA, BUS, BRANCH, GEN, FULLJAC)
0007 %   [J, YBUS, YF, YT] = MAKEJAC(MPC)
0008 %
0009 %   Returns the power flow Jacobian and, optionally, the system admittance
0010 %   matrices. Inputs can be a MATPOWER case struct or individual BASEMVA,
0011 %   BUS, BRANCH and GEN values. Bus numbers must be consecutive beginning
0012 %   at 1 (i.e. internal ordering). If the FULLJAC argument is present and
0013 %   true, it returns the full Jacobian (sensitivities of all bus injections
0014 %   w.r.t all voltage angles/magnitudes) as opposed to the reduced version
0015 %   used in the Newton power flow updates. The units for all quantities are
0016 %   in per unit with radians for voltage angles.
0017 %
0018 %   Note: This function builds the Jacobian from scratch, rebuilding the
0019 %         YBUS matrix in the process. You probably don't want to use this
0020 %         in performance critical code.
0021 %
0022 %   See also MAKEYBUS, EXT2INT
0023 
0024 %   MATPOWER
0025 %   Copyright (c) 1996-2017, Power Systems Engineering Research Center (PSERC)
0026 %   by Ray Zimmerman, PSERC Cornell
0027 %
0028 %   This file is part of MATPOWER.
0029 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0030 %   See https://matpower.org for more info.
0031 
0032 if nargin < 4
0033     mpc     = baseMVA;
0034     if nargin > 1
0035         fullJac = bus;
0036     else
0037         fullJac = 0;
0038     end
0039     baseMVA = mpc.baseMVA;
0040     bus     = mpc.bus;
0041     branch  = mpc.branch;
0042     gen     = mpc.gen;
0043 elseif nargin < 5
0044     fullJac = 0;
0045 end
0046 
0047 %% define named indices into bus, gen, branch matrices
0048 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0049     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0050 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0051     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0052     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0053 
0054 %% build Ybus
0055 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0056 
0057 %% extract voltage
0058 V = bus(:, VM) .* exp(1j * pi/180 * bus(:, VA));
0059 
0060 %% make sure we use generator setpoint voltage for PV and slack buses
0061 on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
0062 gbus = gen(on, GEN_BUS);                %% what buses are they at?
0063 k = find(bus(gbus, BUS_TYPE) == PV | bus(gbus, BUS_TYPE) == REF);
0064 V(gbus(k)) = gen(on(k), VG) ./ abs(V(gbus(k))).* V(gbus(k));
0065 
0066 %% build Jacobian
0067 [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V);
0068 if fullJac
0069     j11 = real(dSbus_dVa);
0070     j12 = real(dSbus_dVm);
0071     j21 = imag(dSbus_dVa);
0072     j22 = imag(dSbus_dVm);
0073 else
0074     %% get bus index lists of each type of bus
0075     [ref, pv, pq] = bustypes(bus, gen);
0076 
0077     j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
0078     j12 = real(dSbus_dVm([pv; pq], pq));
0079     j21 = imag(dSbus_dVa(pq, [pv; pq]));
0080     j22 = imag(dSbus_dVm(pq, pq));
0081 end
0082 
0083 J = [   j11 j12;
0084         j21 j22;    ];