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 % $Id: opf_consfcn.m 2266 2014-01-16 18:12:59Z ray $ 0040 % by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales 0041 % and Ray Zimmerman, PSERC Cornell 0042 % Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC) 0043 % 0044 % This file is part of MATPOWER. 0045 % See http://www.pserc.cornell.edu/matpower/ for more info. 0046 % 0047 % MATPOWER is free software: you can redistribute it and/or modify 0048 % it under the terms of the GNU General Public License as published 0049 % by the Free Software Foundation, either version 3 of the License, 0050 % or (at your option) any later version. 0051 % 0052 % MATPOWER is distributed in the hope that it will be useful, 0053 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0054 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0055 % GNU General Public License for more details. 0056 % 0057 % You should have received a copy of the GNU General Public License 0058 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0059 % 0060 % Additional permission under GNU GPL version 3 section 7 0061 % 0062 % If you modify MATPOWER, or any covered work, to interface with 0063 % other modules (such as MATLAB code and MEX-files) available in a 0064 % MATLAB(R) or comparable environment containing parts covered 0065 % under other licensing terms, the licensors of MATPOWER grant 0066 % you additional permission to convey the resulting work. 0067 0068 %%----- initialize ----- 0069 %% define named indices into data matrices 0070 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0071 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0072 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0073 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0074 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0075 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0076 0077 %% unpack data 0078 mpc = get_mpc(om); 0079 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); 0080 vv = get_idx(om); 0081 0082 %% problem dimensions 0083 nb = size(bus, 1); %% number of buses 0084 nl = size(branch, 1); %% number of branches 0085 ng = size(gen, 1); %% number of dispatchable injections 0086 nxyz = length(x); %% total number of control vars of all types 0087 0088 %% set default constrained lines 0089 if nargin < 7 0090 il = (1:nl); %% all lines have limits by default 0091 end 0092 nl2 = length(il); %% number of constrained lines 0093 0094 %% grab Pg & Qg 0095 Pg = x(vv.i1.Pg:vv.iN.Pg); %% active generation in p.u. 0096 Qg = x(vv.i1.Qg:vv.iN.Qg); %% reactive generation in p.u. 0097 0098 %% put Pg & Qg back in gen 0099 gen(:, PG) = Pg * baseMVA; %% active generation in MW 0100 gen(:, QG) = Qg * baseMVA; %% reactive generation in MVAr 0101 0102 %% rebuild Sbus 0103 Sbus = makeSbus(baseMVA, bus, gen); %% net injected power in p.u. 0104 0105 %% ----- evaluate constraints ----- 0106 %% reconstruct V 0107 Va = x(vv.i1.Va:vv.iN.Va); 0108 Vm = x(vv.i1.Vm:vv.iN.Vm); 0109 V = Vm .* exp(1j * Va); 0110 0111 %% evaluate power flow equations 0112 mis = V .* conj(Ybus * V) - Sbus; 0113 0114 %%----- evaluate constraint function values ----- 0115 %% first, the equality constraints (power flow) 0116 g = [ real(mis); %% active power mismatch for all buses 0117 imag(mis) ]; %% reactive power mismatch for all buses 0118 0119 %% then, the inequality constraints (branch flow limits) 0120 if nl2 > 0 0121 flow_max = (branch(il, RATE_A)/baseMVA).^2; 0122 flow_max(flow_max == 0) = Inf; 0123 if upper(mpopt.opf.flow_lim(1)) == 'I' %% current magnitude limit, |I| 0124 If = Yf * V; 0125 It = Yt * V; 0126 h = [ If .* conj(If) - flow_max; %% branch current limits (from bus) 0127 It .* conj(It) - flow_max ]; %% branch current limits (to bus) 0128 else 0129 %% compute branch power flows 0130 Sf = V(branch(il, F_BUS)) .* conj(Yf * V); %% complex power injected at "from" bus (p.u.) 0131 St = V(branch(il, T_BUS)) .* conj(Yt * V); %% complex power injected at "to" bus (p.u.) 0132 if upper(mpopt.opf.flow_lim(1)) == 'P' %% active power limit, P (Pan Wei) 0133 h = [ real(Sf).^2 - flow_max; %% branch real power limits (from bus) 0134 real(St).^2 - flow_max ]; %% branch real power limits (to bus) 0135 else %% apparent power limit, |S| 0136 h = [ Sf .* conj(Sf) - flow_max; %% branch apparent power limits (from bus) 0137 St .* conj(St) - flow_max ]; %% branch apparent power limits (to bus) 0138 end 0139 end 0140 else 0141 h = zeros(0,1); 0142 end 0143 0144 %%----- evaluate partials of constraints ----- 0145 if nargout > 2 0146 %% index ranges 0147 iVa = vv.i1.Va:vv.iN.Va; 0148 iVm = vv.i1.Vm:vv.iN.Vm; 0149 iPg = vv.i1.Pg:vv.iN.Pg; 0150 iQg = vv.i1.Qg:vv.iN.Qg; 0151 0152 %% compute partials of injected bus powers 0153 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); %% w.r.t. V 0154 neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng); %% Pbus w.r.t. Pg 0155 %% Qbus w.r.t. Qg 0156 0157 %% construct Jacobian of equality (power flow) constraints and transpose it 0158 dg = sparse(2*nb, nxyz); 0159 dg(:, [iVa iVm iPg iQg]) = [ 0160 real([dSbus_dVa dSbus_dVm]) neg_Cg sparse(nb, ng); %% P mismatch w.r.t Va, Vm, Pg, Qg 0161 imag([dSbus_dVa dSbus_dVm]) sparse(nb, ng) neg_Cg; %% Q mismatch w.r.t Va, Vm, Pg, Qg 0162 ]; 0163 dg = dg'; 0164 0165 if nl2 > 0 0166 %% compute partials of Flows w.r.t. V 0167 if upper(mpopt.opf.flow_lim(1)) == 'I' %% current 0168 [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dIbr_dV(branch(il,:), Yf, Yt, V); 0169 else %% power 0170 [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dSbr_dV(branch(il,:), Yf, Yt, V); 0171 end 0172 if upper(mpopt.opf.flow_lim(1)) == 'P' %% real part of flow (active power) 0173 dFf_dVa = real(dFf_dVa); 0174 dFf_dVm = real(dFf_dVm); 0175 dFt_dVa = real(dFt_dVa); 0176 dFt_dVm = real(dFt_dVm); 0177 Ff = real(Ff); 0178 Ft = real(Ft); 0179 end 0180 0181 %% squared magnitude of flow (of complex power or current, or real power) 0182 [df_dVa, df_dVm, dt_dVa, dt_dVm] = ... 0183 dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft); 0184 0185 %% construct Jacobian of inequality (branch flow) constraints & transpose 0186 dh = sparse(2*nl2, nxyz); 0187 dh(:, [iVa iVm]) = [ 0188 df_dVa, df_dVm; %% "from" flow limit 0189 dt_dVa, dt_dVm; %% "to" flow limit 0190 ]; 0191 dh = dh'; 0192 else 0193 dh = sparse(nxyz, 0); 0194 end 0195 end