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];