opf_current_balance_hess

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

   Hessian evaluation function for AC real and imaginary current balance
   constraints.

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

   Outputs:
     D2G : Hessian of current balance constraints.

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

   See also OPF_CURRENT_BALANCE_FCN.

0001 function d2G = opf_current_balance_hess(x, lambda, mpc, Ybus, mpopt)
0002 %OPF_CURRENT_BALANCE_HESS  Evaluates Hessian of current balance constraints.
0003 %   D2G = OPF_CURRENT_BALANCE_HESS(X, LAMBDA, OM, YBUS, MPOPT)
0004 %
0005 %   Hessian evaluation function for AC real and imaginary current balance
0006 %   constraints.
0007 %
0008 %   Inputs:
0009 %     X : optimization vector
0010 %     LAMBDA : column vector of Lagrange multipliers on real and imaginary
0011 %              current balance constraints
0012 %     MPC : MATPOWER case struct
0013 %     YBUS : bus admittance matrix
0014 %     MPOPT : MATPOWER options struct
0015 %
0016 %   Outputs:
0017 %     D2G : Hessian of current balance constraints.
0018 %
0019 %   Example:
0020 %       d2G = opf_current_balance_hess(x, lambda, mpc, Ybus, mpopt);
0021 %
0022 %   See also OPF_CURRENT_BALANCE_FCN.
0023 
0024 %   MATPOWER
0025 %   Copyright (c) 1996-2018, Power Systems Engineering Research Center (PSERC)
0026 %   by Baljinnyam Sereeter, Delft University of Technology
0027 %   and Ray Zimmerman, PSERC Cornell
0028 %
0029 %   This file is part of MATPOWER.
0030 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0031 %   See https://matpower.org for more info.
0032 
0033 %%----- initialize -----
0034 %% define named indices into data matrices
0035 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0036     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0037     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0038 
0039 %% unpack data
0040 [baseMVA, bus, gen] = deal(mpc.baseMVA, mpc.bus, mpc.gen);
0041 if mpopt.opf.v_cartesian
0042     [Vr, Vi, Pg, Qg] = deal(x{:});
0043     V = Vr + 1j * Vi;           %% reconstruct V
0044 else
0045     [Va, Vm, Pg, Qg] = deal(x{:});
0046     V = Vm .* exp(1j * Va);     %% reconstruct V
0047 end
0048 
0049 %% problem dimensions
0050 nb = length(V);             %% number of buses
0051 ng = length(Pg);            %% number of dispatchable injections
0052 
0053 nlam = length(lambda) / 2;
0054 lamP = lambda(1:nlam);
0055 lamQ = lambda((1:nlam)+nlam);
0056 Cg = sparse(gen(:, GEN_BUS), 1:ng, 1, nb, ng);
0057 
0058 %% put Pg, Qg back in gen
0059 gen(:, PG) = Pg * baseMVA;  %% active generation in MW
0060 gen(:, QG) = Qg * baseMVA;  %% reactive generation in MVAr
0061 
0062 %% rebuild Sbus
0063 Sbus = makeSbus(baseMVA, bus, gen);
0064 
0065 %%----- evaluate Hessian of current balance constraints -----
0066 %% compute 2nd derivatives
0067 [Gr11, Gr12, Gr21, Gr22] = d2Imis_dV2(Sbus, Ybus, V, lamP, mpopt.opf.v_cartesian);
0068 [Gi11, Gi12, Gi21, Gi22] = d2Imis_dV2(Sbus, Ybus, V, lamQ, mpopt.opf.v_cartesian);
0069 Gr_sv = d2Imis_dVdSg(Cg, V, lamP, mpopt.opf.v_cartesian);
0070 Gi_sv = d2Imis_dVdSg(Cg, V, lamQ, mpopt.opf.v_cartesian);
0071 
0072 %% construct Hessian
0073 d2G = [
0074     real([Gr11 Gr12; Gr21 Gr22]) + imag([Gi11 Gi12; Gi21 Gi22]), real(Gr_sv.') + imag(Gi_sv.');
0075     real(Gr_sv) + imag(Gi_sv), sparse(2*ng, 2*ng)
0076 ];