OPF_CONSFCN Evaluates nonlinear constraints and their Jacobian for OPF. [H, G, DH, DG] = OPF_CONSFCN(X, OM, YBUS, YF, YT, MPOPT, IL) Constraint evaluation function for AC optimal power flow, suitable for use with MIPS or FMINCON. Computes constraint vectors and their gradients. Inputs: X : optimization vector OM : OPF model object YBUS : bus admittance matrix YF : admittance matrix for "from" end of constrained branches YT : admittance matrix for "to" end of constrained branches MPOPT : MATPOWER options struct IL : (optional) vector of branch indices corresponding to branches with flow limits (all others are assumed to be unconstrained). The default is [1:nl] (all branches). YF and YT contain only the rows corresponding to IL. Outputs: H : vector of inequality constraint values (flow limits) limit^2 - flow^2, where the flow can be apparent power real power or current, depending on value of opf.flow_lim in MPOPT (only for constrained lines) G : vector of equality constraint values (power balances) DH : (optional) inequality constraint gradients, column j is gradient of H(j) DG : (optional) equality constraint gradients Examples: [h, g] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt); [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt); [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt, il); See also OPF_COSTFCN, OPF_HESSFCN.
0001 function [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt, il, varargin) 0002 %OPF_CONSFCN Evaluates nonlinear constraints and their Jacobian for OPF. 0003 % [H, G, DH, DG] = OPF_CONSFCN(X, OM, YBUS, YF, YT, MPOPT, IL) 0004 % 0005 % Constraint evaluation function for AC optimal power flow, suitable 0006 % for use with MIPS or FMINCON. Computes constraint vectors and their 0007 % gradients. 0008 % 0009 % Inputs: 0010 % X : optimization vector 0011 % OM : OPF model object 0012 % YBUS : bus admittance matrix 0013 % YF : admittance matrix for "from" end of constrained branches 0014 % YT : admittance matrix for "to" end of constrained branches 0015 % MPOPT : MATPOWER options struct 0016 % IL : (optional) vector of branch indices corresponding to 0017 % branches with flow limits (all others are assumed to be 0018 % unconstrained). The default is [1:nl] (all branches). 0019 % YF and YT contain only the rows corresponding to IL. 0020 % 0021 % Outputs: 0022 % H : vector of inequality constraint values (flow limits) 0023 % limit^2 - flow^2, where the flow can be apparent power 0024 % real power or current, depending on value of 0025 % opf.flow_lim in MPOPT (only for constrained lines) 0026 % G : vector of equality constraint values (power balances) 0027 % DH : (optional) inequality constraint gradients, column j is 0028 % gradient of H(j) 0029 % DG : (optional) equality constraint gradients 0030 % 0031 % Examples: 0032 % [h, g] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt); 0033 % [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt); 0034 % [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt, il); 0035 % 0036 % See also OPF_COSTFCN, OPF_HESSFCN. 0037 0038 % MATPOWER 0039 % Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) 0040 % by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia 0041 % Ray Zimmerman, PSERC Cornell and Shrirang Abhyankar 0042 % 0043 % This file is part of MATPOWER. 0044 % Covered by the 3-clause BSD License (see LICENSE file for details). 0045 % See http://www.pserc.cornell.edu/matpower/ for more info. 0046 0047 %%----- initialize ----- 0048 %% define named indices into data matrices 0049 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0050 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0051 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0052 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0053 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0054 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0055 0056 %% unpack data 0057 mpc = get_mpc(om); 0058 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); 0059 vv = get_idx(om); 0060 0061 %% problem dimensions 0062 nb = size(bus, 1); %% number of buses 0063 nl = size(branch, 1); %% number of branches 0064 ng = size(gen, 1); %% number of dispatchable injections 0065 nxyz = length(x); %% total number of control vars of all types 0066 0067 %% set default constrained lines 0068 if nargin < 7 0069 il = (1:nl); %% all lines have limits by default 0070 end 0071 nl2 = length(il); %% number of constrained lines 0072 0073 %% grab Pg & Qg 0074 Pg = x(vv.i1.Pg:vv.iN.Pg); %% active generation in p.u. 0075 Qg = x(vv.i1.Qg:vv.iN.Qg); %% reactive generation in p.u. 0076 0077 %% put Pg & Qg back in gen 0078 gen(:, PG) = Pg * baseMVA; %% active generation in MW 0079 gen(:, QG) = Qg * baseMVA; %% reactive generation in MVAr 0080 0081 %% ----- evaluate constraints ----- 0082 %% reconstruct V 0083 Va = x(vv.i1.Va:vv.iN.Va); 0084 Vm = x(vv.i1.Vm:vv.iN.Vm); 0085 V = Vm .* exp(1j * Va); 0086 0087 %% rebuild Sbus 0088 Sbus = makeSbus(baseMVA, bus, gen, mpopt, Vm); %% net injected power in p.u. 0089 0090 %% evaluate power flow equations 0091 mis = V .* conj(Ybus * V) - Sbus; 0092 0093 %%----- evaluate constraint function values ----- 0094 %% first, the equality constraints (power flow) 0095 g = [ real(mis); %% active power mismatch for all buses 0096 imag(mis) ]; %% reactive power mismatch for all buses 0097 0098 %% then, the inequality constraints (branch flow limits) 0099 if nl2 > 0 0100 flow_max = (branch(il, RATE_A)/baseMVA).^2; 0101 flow_max(flow_max == 0) = Inf; 0102 if upper(mpopt.opf.flow_lim(1)) == 'I' %% current magnitude limit, |I| 0103 If = Yf * V; 0104 It = Yt * V; 0105 h = [ If .* conj(If) - flow_max; %% branch current limits (from bus) 0106 It .* conj(It) - flow_max ]; %% branch current limits (to bus) 0107 else 0108 %% compute branch power flows 0109 Sf = V(branch(il, F_BUS)) .* conj(Yf * V); %% complex power injected at "from" bus (p.u.) 0110 St = V(branch(il, T_BUS)) .* conj(Yt * V); %% complex power injected at "to" bus (p.u.) 0111 if upper(mpopt.opf.flow_lim(1)) == 'P' %% active power limit, P (Pan Wei) 0112 h = [ real(Sf).^2 - flow_max; %% branch real power limits (from bus) 0113 real(St).^2 - flow_max ]; %% branch real power limits (to bus) 0114 else %% apparent power limit, |S| 0115 h = [ Sf .* conj(Sf) - flow_max; %% branch apparent power limits (from bus) 0116 St .* conj(St) - flow_max ]; %% branch apparent power limits (to bus) 0117 end 0118 end 0119 else 0120 h = zeros(0,1); 0121 end 0122 0123 %%----- evaluate partials of constraints ----- 0124 if nargout > 2 0125 %% index ranges 0126 iVa = vv.i1.Va:vv.iN.Va; 0127 iVm = vv.i1.Vm:vv.iN.Vm; 0128 iPg = vv.i1.Pg:vv.iN.Pg; 0129 iQg = vv.i1.Qg:vv.iN.Qg; 0130 0131 %% compute partials of injected bus powers 0132 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); %% w.r.t. V 0133 [dummy, neg_dSd_dVm] = makeSbus(baseMVA, bus, gen, mpopt, Vm); 0134 dSbus_dVm = dSbus_dVm - neg_dSd_dVm; 0135 neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng); %% Pbus w.r.t. Pg 0136 %% Qbus w.r.t. Qg 0137 0138 %% construct Jacobian of equality (power flow) constraints and transpose it 0139 dg = sparse(2*nb, nxyz); 0140 dg(:, [iVa iVm iPg iQg]) = [ 0141 real([dSbus_dVa dSbus_dVm]) neg_Cg sparse(nb, ng); %% P mismatch w.r.t Va, Vm, Pg, Qg 0142 imag([dSbus_dVa dSbus_dVm]) sparse(nb, ng) neg_Cg; %% Q mismatch w.r.t Va, Vm, Pg, Qg 0143 ]; 0144 dg = dg'; 0145 0146 if nl2 > 0 0147 %% compute partials of Flows w.r.t. V 0148 if upper(mpopt.opf.flow_lim(1)) == 'I' %% current 0149 [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dIbr_dV(branch(il,:), Yf, Yt, V); 0150 else %% power 0151 [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dSbr_dV(branch(il,:), Yf, Yt, V); 0152 end 0153 if upper(mpopt.opf.flow_lim(1)) == 'P' %% real part of flow (active power) 0154 dFf_dVa = real(dFf_dVa); 0155 dFf_dVm = real(dFf_dVm); 0156 dFt_dVa = real(dFt_dVa); 0157 dFt_dVm = real(dFt_dVm); 0158 Ff = real(Ff); 0159 Ft = real(Ft); 0160 end 0161 0162 %% squared magnitude of flow (of complex power or current, or real power) 0163 [df_dVa, df_dVm, dt_dVa, dt_dVm] = ... 0164 dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft); 0165 0166 %% construct Jacobian of inequality (branch flow) constraints & transpose 0167 dh = sparse(2*nl2, nxyz); 0168 dh(:, [iVa iVm]) = [ 0169 df_dVa, df_dVm; %% "from" flow limit 0170 dt_dVa, dt_dVm; %% "to" flow limit 0171 ]; 0172 dh = dh'; 0173 else 0174 dh = sparse(nxyz, 0); 0175 end 0176 end