OPF_BRANCH_FLOW_FCN Evaluates AC branch flow constraints and Jacobian.
[H, DH] = OPF_BRANCH_FLOW_FCN(X, OM, YF, YT, IL, MPOPT)
Active power balance equality constraints for AC optimal power flow.
Computes constraint vectors and their gradients.
Inputs:
X : optimization vector
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:
H : vector of inequality constraint values (flow limits)
where the flow can be apparent power, real power, or
current, depending on the value of opf.flow_lim in MPOPT
(only for constrained lines), normally expressed as
(limit^2 - flow^2), except when opf.flow_lim == 'P',
in which case it is simply (limit - flow).
DH : (optional) inequality constraint gradients, column j is
gradient of H(j)
Examples:
h = opf_branch_flow_fcn(x, mpc, Yf, Yt, il, mpopt);
[h, dh] = opf_branch_flow_fcn(x, mpc, Yf, Yt, il, mpopt);
See also OPF_BRANCH_FLOW_HESS.0001 function [h, dh] = opf_branch_flow_fcn(x, mpc, Yf, Yt, il, mpopt) 0002 %OPF_BRANCH_FLOW_FCN Evaluates AC branch flow constraints and Jacobian. 0003 % [H, DH] = OPF_BRANCH_FLOW_FCN(X, OM, YF, YT, IL, MPOPT) 0004 % 0005 % Active power balance equality constraints for AC optimal power flow. 0006 % Computes constraint vectors and their gradients. 0007 % 0008 % Inputs: 0009 % X : optimization vector 0010 % MPC : MATPOWER case struct 0011 % YF : admittance matrix for "from" end of constrained branches 0012 % YT : admittance matrix for "to" end of constrained branches 0013 % IL : vector of branch indices corresponding to branches with 0014 % flow limits (all others are assumed to be unconstrained). 0015 % YF and YT contain only the rows corresponding to IL. 0016 % MPOPT : MATPOWER options struct 0017 % 0018 % Outputs: 0019 % H : vector of inequality constraint values (flow limits) 0020 % where the flow can be apparent power, real power, or 0021 % current, depending on the value of opf.flow_lim in MPOPT 0022 % (only for constrained lines), normally expressed as 0023 % (limit^2 - flow^2), except when opf.flow_lim == 'P', 0024 % in which case it is simply (limit - flow). 0025 % DH : (optional) inequality constraint gradients, column j is 0026 % gradient of H(j) 0027 % 0028 % Examples: 0029 % h = opf_branch_flow_fcn(x, mpc, Yf, Yt, il, mpopt); 0030 % [h, dh] = opf_branch_flow_fcn(x, mpc, Yf, Yt, il, mpopt); 0031 % 0032 % See also OPF_BRANCH_FLOW_HESS. 0033 0034 % MATPOWER 0035 % Copyright (c) 1996-2018, Power Systems Engineering Research Center (PSERC) 0036 % by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia 0037 % and Ray Zimmerman, PSERC Cornell 0038 % 0039 % This file is part of MATPOWER. 0040 % Covered by the 3-clause BSD License (see LICENSE file for details). 0041 % See https://matpower.org for more info. 0042 0043 %%----- initialize ----- 0044 %% define named indices into data matrices 0045 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0046 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0047 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0048 0049 %% unpack data 0050 lim_type = upper(mpopt.opf.flow_lim(1)); 0051 branch = mpc.branch; 0052 if mpopt.opf.v_cartesian 0053 [Vr, Vi] = deal(x{:}); 0054 V = Vr + 1j * Vi; %% reconstruct V 0055 else 0056 [Va, Vm] = deal(x{:}); 0057 V = Vm .* exp(1j * Va); %% reconstruct V 0058 end 0059 0060 %% problem dimensions 0061 nb = length(V); %% number of buses 0062 nl2 = length(il); %% number of constrained lines 0063 0064 %% ----- evaluate constraints ----- 0065 if nl2 > 0 0066 flow_max = branch(il, RATE_A) / mpc.baseMVA; 0067 if lim_type ~= 'P' %% typically use square of flow 0068 flow_max = flow_max.^2; 0069 end 0070 if lim_type == 'I' %% current magnitude limit, |I| 0071 If = Yf * V; 0072 It = Yt * V; 0073 h = [ If .* conj(If) - flow_max; %% branch current limits (from bus) 0074 It .* conj(It) - flow_max ]; %% branch current limits (to bus) 0075 else 0076 %% compute branch power flows 0077 Sf = V(branch(il, F_BUS)) .* conj(Yf * V); %% complex power injected at "from" bus (p.u.) 0078 St = V(branch(il, T_BUS)) .* conj(Yt * V); %% complex power injected at "to" bus (p.u.) 0079 if lim_type == '2' %% active power limit, P squared (Pan Wei) 0080 h = [ real(Sf).^2 - flow_max; %% branch real power limits (from bus) 0081 real(St).^2 - flow_max ]; %% branch real power limits (to bus) 0082 elseif lim_type == 'P' %% active power limit, P 0083 h = [ real(Sf) - flow_max; %% branch real power limits (from bus) 0084 real(St) - flow_max ]; %% branch real power limits (to bus 0085 else %% apparent power limit, |S| 0086 h = [ Sf .* conj(Sf) - flow_max; %% branch apparent power limits (from bus) 0087 St .* conj(St) - flow_max ]; %% branch apparent power limits (to bus) 0088 end 0089 end 0090 else 0091 h = zeros(0,1); 0092 end 0093 0094 %%----- evaluate partials of constraints ----- 0095 if nargout > 1 0096 if nl2 > 0 0097 %% compute partials of Flows w.r.t. V 0098 if lim_type == 'I' %% current 0099 [dFf_dV1, dFf_dV2, dFt_dV1, dFt_dV2, Ff, Ft] = dIbr_dV(branch(il,:), Yf, Yt, V, mpopt.opf.v_cartesian); 0100 else %% power 0101 [dFf_dV1, dFf_dV2, dFt_dV1, dFt_dV2, Ff, Ft] = dSbr_dV(branch(il,:), Yf, Yt, V, mpopt.opf.v_cartesian); 0102 end 0103 if lim_type == 'P' || lim_type == '2' %% real part of flow (active power) 0104 dFf_dV1 = real(dFf_dV1); 0105 dFf_dV2 = real(dFf_dV2); 0106 dFt_dV1 = real(dFt_dV1); 0107 dFt_dV2 = real(dFt_dV2); 0108 Ff = real(Ff); 0109 Ft = real(Ft); 0110 end 0111 0112 if lim_type == 'P' 0113 %% active power 0114 [df_dV1, df_dV2, dt_dV1, dt_dV2] = deal(dFf_dV1, dFf_dV2, dFt_dV1, dFt_dV2); 0115 else 0116 %% squared magnitude of flow (of complex power or current, or real power) 0117 [df_dV1, df_dV2, dt_dV1, dt_dV2] = ... 0118 dAbr_dV(dFf_dV1, dFf_dV2, dFt_dV1, dFt_dV2, Ff, Ft); 0119 end 0120 0121 %% construct Jacobian of "from" and "to" branch flow ineq constraints 0122 dh = [ df_dV1 df_dV2; %% "from" flow limit 0123 dt_dV1 dt_dV2 ]; %% "to" flow limit 0124 else 0125 dh = sparse(0, 2*nb); 0126 end 0127 end