Home > matpower5.0 > extras > sdp_pf > insolvablepf_limitQ.m

insolvablepf_limitQ

PURPOSE ^

INSOLVABLEPF_LIMITQ A sufficient condition for power flow insolvability

SYNOPSIS ^

function [insolvable,eta,mineigratio] = insolvablepf_limitQ(mpc,mpopt)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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...

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005