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-2015 by Power System Engineering Research Center (PSERC) 0040 % by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales 0041 % and Ray Zimmerman, PSERC Cornell 0042 % 0043 % $Id: opf_consfcn.m 2644 2015-03-11 19:34:22Z ray $ 0044 % 0045 % This file is part of MATPOWER. 0046 % Covered by the 3-clause BSD License (see LICENSE file for details). 0047 % See http://www.pserc.cornell.edu/matpower/ for more info. 0048 0049 %%----- initialize ----- 0050 %% define named indices into data matrices 0051 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0052 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0053 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0054 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0055 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0056 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0057 0058 %% unpack data 0059 mpc = get_mpc(om); 0060 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); 0061 vv = get_idx(om); 0062 0063 %% problem dimensions 0064 nb = size(bus, 1); %% number of buses 0065 nl = size(branch, 1); %% number of branches 0066 ng = size(gen, 1); %% number of dispatchable injections 0067 nxyz = length(x); %% total number of control vars of all types 0068 0069 %% set default constrained lines 0070 if nargin < 7 0071 il = (1:nl); %% all lines have limits by default 0072 end 0073 nl2 = length(il); %% number of constrained lines 0074 0075 %% grab Pg & Qg 0076 Pg = x(vv.i1.Pg:vv.iN.Pg); %% active generation in p.u. 0077 Qg = x(vv.i1.Qg:vv.iN.Qg); %% reactive generation in p.u. 0078 0079 %% put Pg & Qg back in gen 0080 gen(:, PG) = Pg * baseMVA; %% active generation in MW 0081 gen(:, QG) = Qg * baseMVA; %% reactive generation in MVAr 0082 0083 %% rebuild Sbus 0084 Sbus = makeSbus(baseMVA, bus, gen); %% net injected power in p.u. 0085 0086 %% ----- evaluate constraints ----- 0087 %% reconstruct V 0088 Va = x(vv.i1.Va:vv.iN.Va); 0089 Vm = x(vv.i1.Vm:vv.iN.Vm); 0090 V = Vm .* exp(1j * Va); 0091 0092 %% evaluate power flow equations 0093 mis = V .* conj(Ybus * V) - Sbus; 0094 0095 %%----- evaluate constraint function values ----- 0096 %% first, the equality constraints (power flow) 0097 g = [ real(mis); %% active power mismatch for all buses 0098 imag(mis) ]; %% reactive power mismatch for all buses 0099 0100 %% then, the inequality constraints (branch flow limits) 0101 if nl2 > 0 0102 flow_max = (branch(il, RATE_A)/baseMVA).^2; 0103 flow_max(flow_max == 0) = Inf; 0104 if upper(mpopt.opf.flow_lim(1)) == 'I' %% current magnitude limit, |I| 0105 If = Yf * V; 0106 It = Yt * V; 0107 h = [ If .* conj(If) - flow_max; %% branch current limits (from bus) 0108 It .* conj(It) - flow_max ]; %% branch current limits (to bus) 0109 else 0110 %% compute branch power flows 0111 Sf = V(branch(il, F_BUS)) .* conj(Yf * V); %% complex power injected at "from" bus (p.u.) 0112 St = V(branch(il, T_BUS)) .* conj(Yt * V); %% complex power injected at "to" bus (p.u.) 0113 if upper(mpopt.opf.flow_lim(1)) == 'P' %% active power limit, P (Pan Wei) 0114 h = [ real(Sf).^2 - flow_max; %% branch real power limits (from bus) 0115 real(St).^2 - flow_max ]; %% branch real power limits (to bus) 0116 else %% apparent power limit, |S| 0117 h = [ Sf .* conj(Sf) - flow_max; %% branch apparent power limits (from bus) 0118 St .* conj(St) - flow_max ]; %% branch apparent power limits (to bus) 0119 end 0120 end 0121 else 0122 h = zeros(0,1); 0123 end 0124 0125 %%----- evaluate partials of constraints ----- 0126 if nargout > 2 0127 %% index ranges 0128 iVa = vv.i1.Va:vv.iN.Va; 0129 iVm = vv.i1.Vm:vv.iN.Vm; 0130 iPg = vv.i1.Pg:vv.iN.Pg; 0131 iQg = vv.i1.Qg:vv.iN.Qg; 0132 0133 %% compute partials of injected bus powers 0134 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); %% w.r.t. V 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