opf_branch_flow_hess

OPF_BRANCH_FLOW_HESS  Evaluates Hessian of branch flow constraints.
   D2H = OPF_BRANCH_FLOW_HESS(X, LAMBDA, OM, YF, YT, IL, MPOPT)

   Hessian evaluation function for AC branch flow constraints.

   Inputs:
     X : optimization vector
     LAMBDA : column vector of Kuhn-Tucker multipliers on constrained
              branch flows
     MPC : MATPOWER case struct
     YF : admittance matrix for "from" end of constrained branches
     YT : admittance matrix for "to" end of constrained branches
     IL : vector of branch indices corresponding to branches with
          flow limits (all others are assumed to be unconstrained).
          YF and YT contain only the rows corresponding to IL.
     MPOPT : MATPOWER options struct

   Outputs:
     D2H : Hessian of AC branch flow constraints.

   Example:
       d2H = opf_branch_flow_hess(x, lambda, mpc, Yf, Yt, il, mpopt);

   See also OPF_BRANCH_FLOW_FCN.

0001 function d2H = opf_branch_flow_hess(x, lambda, mpc, Yf, Yt, il, mpopt)
0002 %OPF_BRANCH_FLOW_HESS  Evaluates Hessian of branch flow constraints.
0003 %   D2H = OPF_BRANCH_FLOW_HESS(X, LAMBDA, OM, YF, YT, IL, MPOPT)
0004 %
0005 %   Hessian evaluation function for AC branch flow constraints.
0006 %
0007 %   Inputs:
0008 %     X : optimization vector
0009 %     LAMBDA : column vector of Kuhn-Tucker multipliers on constrained
0010 %              branch flows
0011 %     MPC : MATPOWER case struct
0012 %     YF : admittance matrix for "from" end of constrained branches
0013 %     YT : admittance matrix for "to" end of constrained branches
0014 %     IL : vector of branch indices corresponding to branches with
0015 %          flow limits (all others are assumed to be unconstrained).
0016 %          YF and YT contain only the rows corresponding to IL.
0017 %     MPOPT : MATPOWER options struct
0018 %
0019 %   Outputs:
0020 %     D2H : Hessian of AC branch flow constraints.
0021 %
0022 %   Example:
0023 %       d2H = opf_branch_flow_hess(x, lambda, mpc, Yf, Yt, il, mpopt);
0024 %
0025 %   See also OPF_BRANCH_FLOW_FCN.
0026 
0027 %   MATPOWER
0028 %   Copyright (c) 1996-2018, Power Systems Engineering Research Center (PSERC)
0029 %   by Ray Zimmerman, PSERC Cornell
0030 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0031 %
0032 %   This file is part of MATPOWER.
0033 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0034 %   See https://matpower.org for more info.
0035 
0036 %%----- initialize -----
0037 %% define named indices into data matrices
0038 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0039     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0040     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0041 
0042 %% unpack data
0043 lim_type = upper(mpopt.opf.flow_lim(1));
0044 if mpopt.opf.v_cartesian
0045     [Vr, Vi] = deal(x{:});
0046     V = Vr + 1j * Vi;           %% reconstruct V
0047 else
0048     [Va, Vm] = deal(x{:});
0049     V = Vm .* exp(1j * Va);     %% reconstruct V
0050 end
0051 
0052 %% problem dimensions
0053 nb = length(V);         %% number of buses
0054 nl2 = length(il);       %% number of constrained lines
0055 
0056 %%----- evaluate Hessian of flow constraints -----
0057 %% keep dimensions of empty matrices/vectors compatible
0058 %% (required to avoid problems when using Knitro
0059 %%  on cases with all lines unconstrained)
0060 nmu = length(lambda) / 2;
0061 if nmu
0062     muF = lambda(1:nmu);
0063     muT = lambda((1:nmu)+nmu);
0064 else    %% keep dimensions of empty matrices/vectors compatible
0065     muF = zeros(0,1);   %% (required to avoid problems when using Knitro
0066     muT = zeros(0,1);   %%  on cases with all lines unconstrained)
0067 end
0068 if lim_type == 'I'          %% square of current
0069     [dIf_dV1, dIf_dV2, dIt_dV1, dIt_dV2, If, It] = dIbr_dV(mpc.branch(il,:), Yf, Yt, V, mpopt.opf.v_cartesian);
0070     d2If_dV2 = @(V, mu)d2Ibr_dV2(Yf, V, mu, mpopt.opf.v_cartesian);
0071     d2It_dV2 = @(V, mu)d2Ibr_dV2(Yt, V, mu, mpopt.opf.v_cartesian);
0072     [Hf11, Hf12, Hf21, Hf22] = d2Abr_dV2(d2If_dV2, dIf_dV1, dIf_dV2, If, V, muF);
0073     [Ht11, Ht12, Ht21, Ht22] = d2Abr_dV2(d2It_dV2, dIt_dV1, dIt_dV2, It, V, muT);
0074 else
0075     f = mpc.branch(il, F_BUS);    %% list of "from" buses
0076     t = mpc.branch(il, T_BUS);    %% list of "to" buses
0077     Cf = sparse(1:nl2, f, ones(nl2, 1), nl2, nb);   %% connection matrix for line & from buses
0078     Ct = sparse(1:nl2, t, ones(nl2, 1), nl2, nb);   %% connection matrix for line & to buses
0079     [dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = dSbr_dV(mpc.branch(il,:), Yf, Yt, V, mpopt.opf.v_cartesian);
0080     d2Sf_dV2 = @(V, mu)d2Sbr_dV2(Cf, Yf, V, mu, mpopt.opf.v_cartesian);
0081     d2St_dV2 = @(V, mu)d2Sbr_dV2(Ct, Yt, V, mu, mpopt.opf.v_cartesian);
0082     if lim_type == '2'        %% square of real power
0083         [Hf11, Hf12, Hf21, Hf22] = d2Abr_dV2(d2Sf_dV2, real(dSf_dV1), real(dSf_dV2), real(Sf), V, muF);
0084         [Ht11, Ht12, Ht21, Ht22] = d2Abr_dV2(d2St_dV2, real(dSt_dV1), real(dSt_dV2), real(St), V, muT);
0085     elseif lim_type == 'P'    %% real power
0086         [Hf11, Hf12, Hf21, Hf22] = d2Sf_dV2(V, muF);
0087         [Ht11, Ht12, Ht21, Ht22] = d2St_dV2(V, muT);
0088         [Hf11, Hf12, Hf21, Hf22] = deal(real(Hf11), real(Hf12), real(Hf21), real(Hf22));
0089         [Ht11, Ht12, Ht21, Ht22] = deal(real(Ht11), real(Ht12), real(Ht21), real(Ht22));
0090     else                      %% square of apparent power
0091         [Hf11, Hf12, Hf21, Hf22] = d2Abr_dV2(d2Sf_dV2, dSf_dV1, dSf_dV2, Sf, V, muF);
0092         [Ht11, Ht12, Ht21, Ht22] = d2Abr_dV2(d2St_dV2, dSt_dV1, dSt_dV2, St, V, muT);
0093     end
0094 end
0095 d2H = [Hf11 Hf12; Hf21 Hf22] + [Ht11 Ht12; Ht21 Ht22];