opf_power_balance_hess

OPF_POWER_BALANCE_HESS  Evaluates Hessian of power balance constraints.
   D2G = OPF_POWER_BALANCE_HESS(X, LAMBDA, OM, YBUS, MPOPT)

   Hessian evaluation function for AC active and reactive power balance
   constraints.

   Inputs:
     X : optimization vector
     LAMBDA : column vector of Lagrange multipliers on active and reactive
              power balance constraints
     MPC : MATPOWER case struct
     YBUS : bus admittance matrix
     MPOPT : MATPOWER options struct

   Outputs:
     D2G : Hessian of power balance constraints.

   Example:
       d2G = opf_power_balance_hess(x, lambda, mpc, Ybus, mpopt);

   See also OPF_POWER_BALANCE_FCN.

0001 function d2G = opf_power_balance_hess(x, lambda, mpc, Ybus, mpopt)
0002 %OPF_POWER_BALANCE_HESS  Evaluates Hessian of power balance constraints.
0003 %   D2G = OPF_POWER_BALANCE_HESS(X, LAMBDA, OM, YBUS, MPOPT)
0004 %
0005 %   Hessian evaluation function for AC active and reactive power balance
0006 %   constraints.
0007 %
0008 %   Inputs:
0009 %     X : optimization vector
0010 %     LAMBDA : column vector of Lagrange multipliers on active and reactive
0011 %              power balance constraints
0012 %     MPC : MATPOWER case struct
0013 %     YBUS : bus admittance matrix
0014 %     MPOPT : MATPOWER options struct
0015 %
0016 %   Outputs:
0017 %     D2G : Hessian of power balance constraints.
0018 %
0019 %   Example:
0020 %       d2G = opf_power_balance_hess(x, lambda, mpc, Ybus, mpopt);
0021 %
0022 %   See also OPF_POWER_BALANCE_FCN.
0023 
0024 %   MATPOWER
0025 %   Copyright (c) 1996-2017, Power Systems Engineering Research Center (PSERC)
0026 %   by Ray Zimmerman, PSERC Cornell
0027 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0028 %   and Baljinnyam Sereeter, Delft University of Technology
0029 %
0030 %   This file is part of MATPOWER.
0031 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0032 %   See https://matpower.org for more info.
0033 
0034 %%----- initialize -----
0035 %% unpack data
0036 if mpopt.opf.v_cartesian
0037     [Vr, Vi, Pg, Qg] = deal(x{:});
0038     V = Vr + 1j * Vi;           %% reconstruct V
0039 else
0040     [Va, Vm, Pg, Qg] = deal(x{:});
0041     V = Vm .* exp(1j * Va);     %% reconstruct V
0042 end
0043 
0044 %% problem dimensions
0045 nb = length(V);             %% number of buses
0046 ng = length(Pg);            %% number of dispatchable injections
0047 
0048 nlam = length(lambda) / 2;
0049 lamP = lambda(1:nlam);
0050 lamQ = lambda((1:nlam)+nlam);
0051 
0052 %%----- evaluate Hessian of power balance constraints -----
0053 %% compute 2nd derivatives
0054 [Gp11, Gp12, Gp21, Gp22] = d2Sbus_dV2(Ybus, V, lamP, mpopt.opf.v_cartesian);
0055 [Gq11, Gq12, Gq21, Gq22] = d2Sbus_dV2(Ybus, V, lamQ, mpopt.opf.v_cartesian);
0056 
0057 if ~mpopt.opf.v_cartesian
0058     %% adjust for voltage dependent loads (constant impedance part of ZIP loads)
0059     diaglam = sparse(1:nb, 1:nb, lamP, nb, nb);
0060     Sd = makeSdzip(mpc.baseMVA, mpc.bus, mpopt);
0061     diagSdz = sparse(1:nb, 1:nb, Sd.z, nb, nb);
0062     Gp22 = Gp22 + 2 * diaglam * diagSdz;
0063 end
0064 
0065 %% construct Hessian
0066 d2G = [
0067     real([Gp11 Gp12; Gp21 Gp22]) + imag([Gq11 Gq12; Gq21 Gq22]) sparse(2*nb, 2*ng);
0068     sparse(2*ng, 2*nb + 2*ng)
0069 ];