INSOLVABLEPF_LIMITQ A sufficient condition for power flow insolvability
considering generator reactive power limits
[INSOLVABLE,VSLACK_MIN,SIGMA,ETA,MINEIGRATIO] =
INSOLVABLEPF_LIMITQ(MPC,MPOPT)
Evaluates a sufficient condition for insolvability of the power flow
equations considering generator reactive power limits. This function
uses a mixed-integer semidefinite programming formulation of the power
flow equations with generator reactive power limits that maximizes the
power injections in a uniform, constant power factor profile. eta
indicates the factor by which the power injections can be increased in
this profile. If eta < 1, no power flow solution exists that satisfies
the generator reactive power limits. The converse does not
necessarily hold; a power flow solution may not exist even for cases
that this function does not indicate are insolvable. See [1] for
further details.
Note that this function is only suitable for small systems due to the
computational requirements of the mixed-integer semidefinite
programming solver in YALMIP.
Inputs:
MPC : A MATPOWER case specifying the desired power flow equations.
MPOPT : A MATPOWER options struct. If not specified, it is
assumed to be the default mpoption.
Outputs:
INSOLVABLE : Binary variable. A value of 1 indicates that the
specified power flow equations are insolvable, while a value of
0 means that the insolvability condition is indeterminant (a
solution may or may not exist).
ETA : Power injection margin to the power flow solvability
boundary in the profile of a uniform, constant power factor
change in power injections.
MINEIGRATIO : A measure of satisfaction of the rank relaxation.
Large values indicate satisfaction. (Note that satisfaction of
the rank relaxation is not required for correctness of the
insolvability condition).
[1] D.K. Molzahn, V. Dawar, B.C. Lesieutre, and C.L. DeMarco, "Sufficient
Conditions for Power Flow Insolvability Considering Reactive Power
Limited Generators with Applications to Voltage Stability Margins,"
in Bulk Power System Dynamics and Control - IX. Optimization,
Security and Control of the Emerging Power Grid, 2013 IREP Symposium,
25-30 August 2013.0001 function [insolvable,eta,mineigratio] = insolvablepf_limitQ(mpc,mpopt) 0002 %INSOLVABLEPF_LIMITQ A sufficient condition for power flow insolvability 0003 %considering generator reactive power limits 0004 % 0005 % [INSOLVABLE,VSLACK_MIN,SIGMA,ETA,MINEIGRATIO] = 0006 % INSOLVABLEPF_LIMITQ(MPC,MPOPT) 0007 % 0008 % Evaluates a sufficient condition for insolvability of the power flow 0009 % equations considering generator reactive power limits. This function 0010 % uses a mixed-integer semidefinite programming formulation of the power 0011 % flow equations with generator reactive power limits that maximizes the 0012 % power injections in a uniform, constant power factor profile. eta 0013 % indicates the factor by which the power injections can be increased in 0014 % this profile. If eta < 1, no power flow solution exists that satisfies 0015 % the generator reactive power limits. The converse does not 0016 % necessarily hold; a power flow solution may not exist even for cases 0017 % that this function does not indicate are insolvable. See [1] for 0018 % further details. 0019 % 0020 % Note that this function is only suitable for small systems due to the 0021 % computational requirements of the mixed-integer semidefinite 0022 % programming solver in YALMIP. 0023 % 0024 % Inputs: 0025 % MPC : A MATPOWER case specifying the desired power flow equations. 0026 % MPOPT : A MATPOWER options struct. If not specified, it is 0027 % assumed to be the default mpoption. 0028 % 0029 % Outputs: 0030 % INSOLVABLE : Binary variable. A value of 1 indicates that the 0031 % specified power flow equations are insolvable, while a value of 0032 % 0 means that the insolvability condition is indeterminant (a 0033 % solution may or may not exist). 0034 % ETA : Power injection margin to the power flow solvability 0035 % boundary in the profile of a uniform, constant power factor 0036 % change in power injections. 0037 % MINEIGRATIO : A measure of satisfaction of the rank relaxation. 0038 % Large values indicate satisfaction. (Note that satisfaction of 0039 % the rank relaxation is not required for correctness of the 0040 % insolvability condition). 0041 % 0042 % [1] D.K. Molzahn, V. Dawar, B.C. Lesieutre, and C.L. DeMarco, "Sufficient 0043 % Conditions for Power Flow Insolvability Considering Reactive Power 0044 % Limited Generators with Applications to Voltage Stability Margins," 0045 % in Bulk Power System Dynamics and Control - IX. Optimization, 0046 % Security and Control of the Emerging Power Grid, 2013 IREP Symposium, 0047 % 25-30 August 2013. 0048 0049 % MATPOWER 0050 % $Id: insolvablepf_limitQ.m 2280 2014-01-17 23:28:37Z ray $ 0051 % by Daniel Molzahn, PSERC U of Wisc, Madison 0052 % Copyright (c) 2013-2014 by Power System Engineering Research Center (PSERC) 0053 % 0054 % This file is part of MATPOWER. 0055 % See http://www.pserc.cornell.edu/matpower/ for more info. 0056 % 0057 % MATPOWER is free software: you can redistribute it and/or modify 0058 % it under the terms of the GNU General Public License as published 0059 % by the Free Software Foundation, either version 3 of the License, 0060 % or (at your option) any later version. 0061 % 0062 % MATPOWER is distributed in the hope that it will be useful, 0063 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0064 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0065 % GNU General Public License for more details. 0066 % 0067 % You should have received a copy of the GNU General Public License 0068 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0069 % 0070 % Additional permission under GNU GPL version 3 section 7 0071 % 0072 % If you modify MATPOWER, or any covered work, to interface with 0073 % other modules (such as MATLAB code and MEX-files) available in a 0074 % MATLAB(R) or comparable environment containing parts covered 0075 % under other licensing terms, the licensors of MATPOWER grant 0076 % you additional permission to convey the resulting work. 0077 0078 if nargin < 2 0079 mpopt = mpoption; 0080 end 0081 0082 mpc = loadcase(mpc); 0083 0084 % Unpack options 0085 ignore_angle_lim = mpopt.opf.ignore_angle_lim; 0086 verbose = mpopt.verbose; 0087 ndisplay = mpopt.sdp_pf.ndisplay; %% Determine display frequency of diagonastic information 0088 0089 maxSystemSize = 57; 0090 fixPVbusInjection = 0; % If equal to 1, don't allow changes in active power injections at PV buses. 0091 0092 if ~have_fcn('yalmip') 0093 error('insolvablepf_limitQ: The software package YALMIP is required to run insolvablepf_limitQ. See http://users.isy.liu.se/johanl/yalmip/'); 0094 end 0095 0096 % set YALMIP options struct in SDP_PF (for details, see help sdpsettings) 0097 sdpopts = yalmip_options([], mpopt); 0098 0099 % Change solver to YALMIP's branch-and-bound algorithm 0100 sdpopts = sdpsettings(sdpopts,'solver','bnb','bnb.solver',sdpopts.solver); 0101 0102 if strcmp(sdpopts.solver, 'sedumi') || strcmp(sdpopts.solver,'sdpt3') 0103 sdpopts = sdpsettings(sdpopts,'solver','bnb'); 0104 end 0105 0106 %% define named indices into data matrices 0107 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0108 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0109 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0110 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0111 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0112 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0113 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0114 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0115 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost; 0116 0117 0118 %% Load mpc data 0119 0120 mpc = ext2int(mpc); 0121 0122 if toggle_dcline(mpc, 'status') 0123 error('insolvablepf_limitQ: DC lines are not implemented in insolvablepf'); 0124 end 0125 0126 nbus = size(mpc.bus,1); 0127 ngen = size(mpc.gen,1); 0128 0129 if nbus > maxSystemSize 0130 error('insolvablepf_limitQ: System is too large (more than 57 buses) to solve with insolvablepf_limitQ'); 0131 end 0132 0133 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360)) 0134 warning('insolvablepf_limitQ: Angle difference constraints are not implemented in SDP_PF. Ignoring angle difference constraints.'); 0135 end 0136 0137 % Some of the larger system (e.g., case2746wp) have generators 0138 % corresponding to buses that have bustype == PQ. Change these 0139 % to PV buses. 0140 for i=1:ngen 0141 busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS)); 0142 if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF) 0143 mpc.bus(busidx,BUS_TYPE) = PV; 0144 if verbose >= 1 0145 warning('insolvablepf_limitQ: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx)); 0146 end 0147 end 0148 end 0149 0150 % Buses may be listed as PV buses without associated generators. Change 0151 % these buses to PQ. 0152 for i=1:nbus 0153 if mpc.bus(i,BUS_TYPE) == PV 0154 genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1); 0155 if isempty(genidx) 0156 mpc.bus(i,BUS_TYPE) = PQ; 0157 if verbose >= 1 0158 warning('insolvablepf_limitQ: PV bus %i has no associated generator! Changing these buses to PQ.',i); 0159 end 0160 end 0161 end 0162 end 0163 0164 0165 Sbase = mpc.baseMVA; 0166 0167 % Create vectors of power injections and voltage magnitudes 0168 Qd = mpc.bus(:,QD) / Sbase; 0169 Qinj = -Qd; 0170 Vmag = mpc.bus(:,VM); 0171 0172 Pd = mpc.bus(:,PD) / Sbase; 0173 Pg = zeros(nbus,1); 0174 Qmin = zeros(nbus,1); 0175 Qmax = zeros(nbus,1); 0176 for i=1:nbus 0177 genidx = find(mpc.gen(:,GEN_BUS) == i); 0178 if ~isempty(genidx) 0179 Pg(i) = sum(mpc.gen(genidx,PG)) / Sbase; 0180 Vmag(i) = mpc.gen(genidx(1),VG); 0181 Qmin(i) = sum(mpc.gen(genidx,QMIN)) / Sbase - Qd(i); 0182 Qmax(i) = sum(mpc.gen(genidx,QMAX)) / Sbase - Qd(i); 0183 end 0184 end 0185 Pinj = Pg - Pd; 0186 0187 0188 %% Functions to build matrices 0189 0190 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc); 0191 0192 %% Create primal SDP variables 0193 0194 [junk1,uniqueGenIdx,junk2] = unique(mpc.gen(:,GEN_BUS)); 0195 mpc.gen = mpc.gen(uniqueGenIdx,:); 0196 ngen = size(mpc.gen,1); 0197 0198 W = sdpvar(2*nbus,2*nbus); 0199 eta = sdpvar(1,1); 0200 constraints = []; 0201 0202 % Binary variables 0203 yL = binvar(ngen,1); 0204 yU = binvar(ngen,1); 0205 0206 % We need a number greater than any plausible value of V^2, but too large 0207 % of a value causes numerical problems 0208 bigM = 10*max(Vmag)^2; 0209 0210 %% Build primal problem 0211 0212 for k=1:nbus 0213 0214 if ~fixPVbusInjection 0215 % PQ and PV buses have uniform active power injection changes 0216 if mpc.bus(k,BUS_TYPE) == PQ || mpc.bus(k,BUS_TYPE) == PV 0217 constraints = [constraints; 0218 trace(Yk(k)*W) == eta*Pinj(k)]; 0219 end 0220 else 0221 % Alternatively, we can fix PV bus active power injections for a 0222 % different power injection profile. 0223 if mpc.bus(k,BUS_TYPE) == PQ 0224 constraints = [constraints; 0225 trace(Yk(k)*W) == eta*Pinj(k)]; 0226 end 0227 if mpc.bus(k,BUS_TYPE) == PV 0228 constraints = [constraints; 0229 trace(Yk(k)*W) == Pinj(k)]; 0230 end 0231 end 0232 0233 % PQ buses have uniform reactive power injection changes 0234 if mpc.bus(k,BUS_TYPE) == PQ 0235 constraints = [constraints; 0236 trace(Yk_(k)*W) == eta*Qinj(k)]; 0237 end 0238 0239 % Mixed integer formulation for generator reactive power limits 0240 if mpc.bus(k,BUS_TYPE) == PV || mpc.bus(k,BUS_TYPE) == REF 0241 0242 genidx = find(mpc.gen(:,GEN_BUS) == k,1); 0243 0244 constraints = [constraints; 0245 trace(Yk_(k)*W) <= Qmin(k)*yL(genidx) + Qmax(k)*(1-yL(genidx)); 0246 trace(Yk_(k)*W) >= Qmax(k)*yU(genidx) + Qmin(k)*(1-yU(genidx))]; 0247 0248 constraints = [constraints; 0249 trace(Mk(k)*W) >= Vmag(k)^2*(1-yU(genidx)); 0250 trace(Mk(k)*W) <= Vmag(k)^2*(1-yL(genidx)) + bigM*(yL(genidx))]; 0251 0252 constraints = [constraints; 0253 yU(genidx) + yL(genidx) <= 1]; 0254 0255 end 0256 0257 if verbose >= 2 && mod(k,ndisplay) == 0 0258 fprintf('SDP creation: Bus %i of %i\n',k,nbus); 0259 end 0260 0261 end % Loop through all buses 0262 0263 % Force at least one bus to supply reactive power balance 0264 % genbuses = (mpc.bus(:,BUS_TYPE) == PV | mpc.bus(:,BUS_TYPE) == REF); 0265 % constraints = [constraints; 0266 % sum(yL(genbuses)+yU(genbuses)) <= ngen - 1]; 0267 constraints = [constraints; 0268 sum(yL + yU) <= ngen - 1]; 0269 0270 constraints = [constraints; W >= 0]; 0271 0272 0273 0274 %% Solve the SDP 0275 0276 % Preserve warning settings 0277 S = warning; 0278 0279 % Run sdp solver 0280 sdpinfo = solvesdp(constraints, -eta, sdpopts); 0281 warning(S); 0282 0283 if verbose >= 2 0284 fprintf('Solver exit message: %s\n',sdpinfo.info); 0285 end 0286 0287 0288 %% Calculate rank characteristics of the solution 0289 0290 evl = sort(eig(double(W))); 0291 mineigratio = abs(evl(end) / evl(end-2)); 0292 0293 0294 %% Output results 0295 0296 eta = double(eta); 0297 insolvable = eta < 1; % Is this right? I think it should be the other way around...