insolvablepfsos

INSOLVABLEPFSOS A sufficient condition for power flow insolvability
using sum of squares programming

   [INSOLVABLE] = INSOLVABLEPFSOS(MPC,MPOPT)

   Uses sum of squares programming to generate an infeasibility 
   certificate for the power flow equations. An infeasibility certificate
   proves insolvability of the power flow equations.

   Note that absence of an infeasibility certificate does not necessarily
   mean that a solution does exist (that is, insolvable = 0 does not imply
   solution existence). See [1] for more details.

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

   [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 = insolvablepfsos(mpc,mpopt)
0002 %INSOLVABLEPFSOS A sufficient condition for power flow insolvability
0003 %using sum of squares programming
0004 %
0005 %   [INSOLVABLE] = INSOLVABLEPFSOS(MPC,MPOPT)
0006 %
0007 %   Uses sum of squares programming to generate an infeasibility
0008 %   certificate for the power flow equations. An infeasibility certificate
0009 %   proves insolvability of the power flow equations.
0010 %
0011 %   Note that absence of an infeasibility certificate does not necessarily
0012 %   mean that a solution does exist (that is, insolvable = 0 does not imply
0013 %   solution existence). See [1] for more details.
0014 %
0015 %   Inputs:
0016 %       MPC : A MATPOWER case specifying the desired power flow equations.
0017 %       MPOPT : A MATPOWER options struct. If not specified, it is
0018 %           assumed to be the default mpoption.
0019 %
0020 %   Outputs:
0021 %       INSOLVABLE : Binary variable. A value of 1 indicates that the
0022 %           specified power flow equations are insolvable, while a value of
0023 %           0 means that the insolvability condition is indeterminant (a
0024 %           solution may or may not exist).
0025 %
0026 %   [1] D.K. Molzahn, V. Dawar, B.C. Lesieutre, and C.L. DeMarco, "Sufficient
0027 %       Conditions for Power Flow Insolvability Considering Reactive Power
0028 %       Limited Generators with Applications to Voltage Stability Margins,"
0029 %       in Bulk Power System Dynamics and Control - IX. Optimization,
0030 %       Security and Control of the Emerging Power Grid, 2013 IREP Symposium,
0031 %       25-30 August 2013.
0032 
0033 %   MATPOWER
0034 %   Copyright (c) 2013-2015 by Power System Engineering Research Center (PSERC)
0035 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0036 %
0037 %   $Id: insolvablepfsos.m 2644 2015-03-11 19:34:22Z ray $
0038 %
0039 %   This file is part of MATPOWER.
0040 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0041 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0042 
0043 define_constants;
0044 
0045 if nargin < 2
0046     mpopt = mpoption;
0047 end
0048 
0049 % Unpack options
0050 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0051 verbose             = mpopt.verbose;
0052 enforce_Qlimits     = mpopt.pf.enforce_q_lims;
0053 
0054 if enforce_Qlimits > 0
0055     enforce_Qlimits = 1;
0056 end
0057 
0058 if ~have_fcn('yalmip')
0059     error('insolvablepfsos: The software package YALMIP is required to run insolvablepfsos. See http://users.isy.liu.se/johanl/yalmip/');
0060 end
0061 
0062 % set YALMIP options struct in SDP_PF (for details, see help sdpsettings)
0063 sdpopts = yalmip_options([], mpopt);
0064 
0065 %% Handle generator reactive power limits
0066 % If no generator reactive power limits are specified, use this code
0067 % directly. If generator reactive power limits are to be enforced, use a
0068 % function that considers reactive power limited generators.
0069 
0070 if enforce_Qlimits
0071     if verbose > 0
0072         fprintf('Generator reactive power limits are enforced. Using function insolvablepfsos_limitQ.\n');
0073     end
0074     insolvable = insolvablepfsos_limitQ(mpc,mpopt);
0075     return;
0076 end
0077 
0078 %% Load data
0079 
0080 mpc = loadcase(mpc);
0081 mpc = ext2int(mpc);
0082 
0083 if toggle_dcline(mpc, 'status')
0084     error('insolvablepfsos: DC lines are not implemented in insolvablepf');
0085 end
0086 
0087 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0088     warning('insolvablepfsos: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0089 end
0090 
0091 nbus = size(mpc.bus,1);
0092 ngen = size(mpc.gen,1);
0093 
0094 % Some of the larger system (e.g., case2746wp) have generators
0095 % corresponding to buses that have bustype == PQ. Change these
0096 % to PV buses.
0097 for i=1:ngen
0098     busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0099     if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0100         mpc.bus(busidx,BUS_TYPE) = PV;
0101         if verbose >= 1
0102             warning('insolvablepfsos: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0103         end
0104     end
0105 end
0106 
0107 % Buses may be listed as PV buses without associated generators. Change
0108 % these buses to PQ.
0109 for i=1:nbus
0110     if mpc.bus(i,BUS_TYPE) == PV
0111         genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0112         if isempty(genidx)
0113             mpc.bus(i,BUS_TYPE) = PQ;
0114             if verbose >= 1
0115                 warning('insolvablepfsos: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0116             end
0117         end
0118     end
0119 end
0120 
0121 pq = find(mpc.bus(:,2) == PQ);
0122 pv = find(mpc.bus(:,2) == PV);
0123 ref = find(mpc.bus(:,2) == REF);
0124 
0125 refpv = sort([ref; pv]);
0126 pvpq = sort([pv; pq]);
0127 
0128 Y = makeYbus(mpc);
0129 G = real(Y);
0130 B = imag(Y);
0131 
0132 % Make vectors of power injections and voltage magnitudes
0133 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0134 P = real(S);
0135 Q = imag(S);
0136 
0137 Vmag = zeros(nbus,1);
0138 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0139 
0140 
0141 %% Create voltage variables
0142 % We need a Vd and Vq for each bus
0143 
0144 if verbose >= 2
0145     fprintf('Creating variables.\n');
0146 end
0147 
0148 Vd = sdpvar(nbus,1);
0149 Vq = sdpvar(nbus,1);
0150 
0151 V = [1; Vd; Vq];
0152 
0153 %% Create polynomials
0154 
0155 if verbose >= 2
0156     fprintf('Creating polynomials.\n');
0157 end
0158 
0159 deg = 0;
0160 
0161 % Polynomials for active power
0162 for i=1:nbus-1
0163     [l,lc] = polynomial(V,deg,0);
0164     lambda(i) = l;
0165     clambda{i} = lc;
0166 end
0167 
0168 % Polynomials for reactive power
0169 for i=1:length(pq)
0170     [g,gc] = polynomial(V,deg,0);
0171     gamma(i) = g;
0172     cgamma{i} = gc;
0173 end
0174 
0175 % Polynomials for voltage magnitude
0176 for i=1:length(refpv)
0177     [m,mc] = polynomial(V,deg,0);
0178     mu(i) = m;
0179     cmu{i} = mc;
0180 end
0181 
0182 % Polynomial for reference angle
0183 [delta, cdelta] = polynomial(V,deg,0);
0184 
0185 
0186 %% Create power flow equation polynomials
0187 
0188 if verbose >= 2
0189     fprintf('Creating power flow equations.\n');
0190 end
0191 
0192 % Make P equations
0193 for k=1:nbus-1
0194     i = pvpq(k);
0195     
0196     % Take advantage of sparsity in forming polynomials
0197 %     fp(k) = Vd(i) * sum(G(:,i).*Vd - B(:,i).*Vq) + Vq(i) * sum(B(:,i).*Vd + G(:,i).*Vq);
0198     
0199     y = find(G(:,i) | B(:,i));
0200     for m=1:length(y)
0201         if exist('fp','var') && length(fp) == k
0202             fp(k) = fp(k) + Vd(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m))) + Vq(i) * (B(y(m),i)*Vd(y(m)) + G(y(m),i)*Vq(y(m)));
0203         else
0204             fp(k) = Vd(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m))) + Vq(i) * (B(y(m),i)*Vd(y(m)) + G(y(m),i)*Vq(y(m)));
0205         end
0206     end
0207 
0208 end
0209 
0210 % Make Q equations
0211 for k=1:length(pq)
0212     i = pq(k);
0213     
0214     % Take advantage of sparsity in forming polynomials
0215 %     fq(k) = Vd(i) * sum(-B(:,i).*Vd - G(:,i).*Vq) + Vq(i) * sum(G(:,i).*Vd - B(:,i).*Vq);
0216     
0217     y = find(G(:,i) | B(:,i));
0218     for m=1:length(y)
0219         if exist('fq','var') && length(fq) == k
0220             fq(k) = fq(k) + Vd(i) * (-B(y(m),i)*Vd(y(m)) - G(y(m),i)*Vq(y(m))) + Vq(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m)));
0221         else
0222             fq(k) = Vd(i) * (-B(y(m),i)*Vd(y(m)) - G(y(m),i)*Vq(y(m))) + Vq(i) * (G(y(m),i)*Vd(y(m)) - B(y(m),i)*Vq(y(m)));
0223         end
0224     end
0225 end
0226 
0227 % Make V equations
0228 for k=1:length(refpv)
0229     i = refpv(k);
0230     fv(k) = Vd(i)^2 + Vq(i)^2;
0231 end
0232 
0233 
0234 %% Form the test polynomial
0235 
0236 if verbose >= 2
0237     fprintf('Forming the test polynomial.\n');
0238 end
0239 
0240 sospoly = 1;
0241 
0242 % Add active power terms
0243 for i=1:nbus-1
0244     sospoly = sospoly + lambda(i) * (fp(i) - P(pvpq(i)));
0245 end
0246 
0247 % Add reactive power terms
0248 for i=1:length(pq)
0249     sospoly = sospoly + gamma(i) * (fq(i) - Q(pq(i)));
0250 end
0251 
0252 % Add voltage magnitude terms
0253 for i=1:length(refpv)
0254     sospoly = sospoly + mu(i) * (fv(i) - Vmag(refpv(i))^2);
0255 end
0256 
0257 % Add angle reference term
0258 sospoly = sospoly + delta * Vq(ref);
0259 
0260 sospoly = -sospoly;
0261 
0262 %% Solve the SOS problem
0263 
0264 if verbose >= 2
0265     fprintf('Solving the sum of squares problem.\n');
0266 end
0267 
0268 sosopts = sdpsettings(sdpopts,'sos.newton',1,'sos.congruence',1);
0269 
0270 constraints = sos(sospoly);
0271 
0272 pvec = cdelta(:).';
0273 
0274 for i=1:length(lambda)
0275     pvec = [pvec clambda{i}(:).'];
0276 end
0277 
0278 for i=1:length(gamma)
0279     pvec = [pvec cgamma{i}(:).'];
0280 end
0281 
0282 for i=1:length(mu)
0283     pvec = [pvec cmu{i}(:).'];
0284 end
0285 
0286 % Preserve warning settings
0287 S = warning;
0288 
0289 sol = solvesos(constraints,[],sosopts,pvec);
0290 
0291 warning(S);
0292 
0293 if verbose >= 2
0294     fprintf('Solver exit message: %s\n',sol.info);
0295 end
0296 
0297 if sol.problem
0298     % Power flow equations may have a solution
0299     insolvable = 0;
0300 elseif sol.problem == 0
0301     % Power flow equations do not have a solution.
0302     insolvable = 1;
0303 end