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-2019, Power Systems Engineering Research Center (PSERC)
0035 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0036 %
0037 %   This file is part of MATPOWER/mx-sdp_pf.
0038 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0039 %   See https://github.com/MATPOWER/mx-sdp_pf/ for more info.
0040 
0041 define_constants;
0042 
0043 if nargin < 2
0044     mpopt = mpoption;
0045 end
0046 
0047 % Unpack options
0048 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0049 verbose             = mpopt.verbose;
0050 enforce_Qlimits     = mpopt.pf.enforce_q_lims;
0051 
0052 if enforce_Qlimits > 0
0053     enforce_Qlimits = 1;
0054 end
0055 
0056 if ~have_feature('yalmip')
0057     error('insolvablepfsos: The software package YALMIP is required to run insolvablepfsos. See https://yalmip.github.io');
0058 end
0059 
0060 % set YALMIP options struct in SDP_PF (for details, see help sdpsettings)
0061 sdpopts = yalmip_options([], mpopt);
0062 
0063 %% Handle generator reactive power limits
0064 % If no generator reactive power limits are specified, use this code
0065 % directly. If generator reactive power limits are to be enforced, use a
0066 % function that considers reactive power limited generators.
0067 
0068 if enforce_Qlimits
0069     if verbose > 0
0070         fprintf('Generator reactive power limits are enforced. Using function insolvablepfsos_limitQ.\n');
0071     end
0072     insolvable = insolvablepfsos_limitQ(mpc,mpopt);
0073     return;
0074 end
0075 
0076 %% Load data
0077 
0078 mpc = loadcase(mpc);
0079 mpc = ext2int(mpc);
0080 
0081 if toggle_dcline(mpc, 'status')
0082     error('insolvablepfsos: DC lines are not implemented in insolvablepf');
0083 end
0084 
0085 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0086     warning('insolvablepfsos: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0087 end
0088 
0089 nbus = size(mpc.bus,1);
0090 ngen = size(mpc.gen,1);
0091 
0092 % Some of the larger system (e.g., case2746wp) have generators
0093 % corresponding to buses that have bustype == PQ. Change these
0094 % to PV buses.
0095 for i=1:ngen
0096     busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0097     if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0098         mpc.bus(busidx,BUS_TYPE) = PV;
0099         if verbose >= 1
0100             warning('insolvablepfsos: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0101         end
0102     end
0103 end
0104 
0105 % Buses may be listed as PV buses without associated generators. Change
0106 % these buses to PQ.
0107 for i=1:nbus
0108     if mpc.bus(i,BUS_TYPE) == PV
0109         genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0110         if isempty(genidx)
0111             mpc.bus(i,BUS_TYPE) = PQ;
0112             if verbose >= 1
0113                 warning('insolvablepfsos: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0114             end
0115         end
0116     end
0117 end
0118 
0119 pq = find(mpc.bus(:,2) == PQ);
0120 pv = find(mpc.bus(:,2) == PV);
0121 ref = find(mpc.bus(:,2) == REF);
0122 
0123 refpv = sort([ref; pv]);
0124 pvpq = sort([pv; pq]);
0125 
0126 Y = makeYbus(mpc);
0127 G = real(Y);
0128 B = imag(Y);
0129 
0130 % Make vectors of power injections and voltage magnitudes
0131 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0132 P = real(S);
0133 Q = imag(S);
0134 
0135 Vmag = zeros(nbus,1);
0136 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0137 
0138 
0139 %% Create voltage variables
0140 % We need a Vd and Vq for each bus
0141 
0142 if verbose >= 2
0143     fprintf('Creating variables.\n');
0144 end
0145 
0146 Vd = sdpvar(nbus,1);
0147 Vq = sdpvar(nbus,1);
0148 
0149 V = [1; Vd; Vq];
0150 
0151 %% Create polynomials
0152 
0153 if verbose >= 2
0154     fprintf('Creating polynomials.\n');
0155 end
0156 
0157 deg = 0;
0158 
0159 % Polynomials for active power
0160 for i=1:nbus-1
0161     [l,lc] = polynomial(V,deg,0);
0162     lambda(i) = l;
0163     clambda{i} = lc;
0164 end
0165 
0166 % Polynomials for reactive power
0167 for i=1:length(pq)
0168     [g,gc] = polynomial(V,deg,0);
0169     gamma(i) = g;
0170     cgamma{i} = gc;
0171 end
0172 
0173 % Polynomials for voltage magnitude
0174 for i=1:length(refpv)
0175     [m,mc] = polynomial(V,deg,0);
0176     mu(i) = m;
0177     cmu{i} = mc;
0178 end
0179 
0180 % Polynomial for reference angle
0181 [delta, cdelta] = polynomial(V,deg,0);
0182 
0183 
0184 %% Create power flow equation polynomials
0185 
0186 if verbose >= 2
0187     fprintf('Creating power flow equations.\n');
0188 end
0189 
0190 % Make P equations
0191 for k=1:nbus-1
0192     i = pvpq(k);
0193     
0194     % Take advantage of sparsity in forming polynomials
0195 %     fp(k) = Vd(i) * sum(G(:,i).*Vd - B(:,i).*Vq) + Vq(i) * sum(B(:,i).*Vd + G(:,i).*Vq);
0196     
0197     y = find(G(:,i) | B(:,i));
0198     for m=1:length(y)
0199         if exist('fp','var') && length(fp) == k
0200             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)));
0201         else
0202             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         end
0204     end
0205 
0206 end
0207 
0208 % Make Q equations
0209 for k=1:length(pq)
0210     i = pq(k);
0211     
0212     % Take advantage of sparsity in forming polynomials
0213 %     fq(k) = Vd(i) * sum(-B(:,i).*Vd - G(:,i).*Vq) + Vq(i) * sum(G(:,i).*Vd - B(:,i).*Vq);
0214     
0215     y = find(G(:,i) | B(:,i));
0216     for m=1:length(y)
0217         if exist('fq','var') && length(fq) == k
0218             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)));
0219         else
0220             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         end
0222     end
0223 end
0224 
0225 % Make V equations
0226 for k=1:length(refpv)
0227     i = refpv(k);
0228     fv(k) = Vd(i)^2 + Vq(i)^2;
0229 end
0230 
0231 
0232 %% Form the test polynomial
0233 
0234 if verbose >= 2
0235     fprintf('Forming the test polynomial.\n');
0236 end
0237 
0238 sospoly = 1;
0239 
0240 % Add active power terms
0241 for i=1:nbus-1
0242     sospoly = sospoly + lambda(i) * (fp(i) - P(pvpq(i)));
0243 end
0244 
0245 % Add reactive power terms
0246 for i=1:length(pq)
0247     sospoly = sospoly + gamma(i) * (fq(i) - Q(pq(i)));
0248 end
0249 
0250 % Add voltage magnitude terms
0251 for i=1:length(refpv)
0252     sospoly = sospoly + mu(i) * (fv(i) - Vmag(refpv(i))^2);
0253 end
0254 
0255 % Add angle reference term
0256 sospoly = sospoly + delta * Vq(ref);
0257 
0258 sospoly = -sospoly;
0259 
0260 %% Solve the SOS problem
0261 
0262 if verbose >= 2
0263     fprintf('Solving the sum of squares problem.\n');
0264 end
0265 
0266 sosopts = sdpsettings(sdpopts,'sos.newton',1,'sos.congruence',1);
0267 
0268 constraints = sos(sospoly);
0269 
0270 pvec = cdelta(:).';
0271 
0272 for i=1:length(lambda)
0273     pvec = [pvec clambda{i}(:).'];
0274 end
0275 
0276 for i=1:length(gamma)
0277     pvec = [pvec cgamma{i}(:).'];
0278 end
0279 
0280 for i=1:length(mu)
0281     pvec = [pvec cmu{i}(:).'];
0282 end
0283 
0284 % Preserve warning settings
0285 S = warning;
0286 if have_feature('matlab', 'vnum') >= 8.006 && have_feature('cplex') && ...
0287         have_feature('cplex', 'vnum') <= 12.006003
0288     warning('OFF', 'MATLAB:lang:badlyScopedReturnValue');
0289 end
0290 
0291 sol = solvesos(constraints,[],sosopts,pvec);
0292 
0293 if ~have_feature('octave') || have_feature('octave', 'vnum') >= 4.001
0294     %% (avoid bug in Octave 4.0.x, where warning state is left corrupted)
0295     warning(S);
0296 end
0297 
0298 if verbose >= 2
0299     fprintf('Solver exit message: %s\n',sol.info);
0300 end
0301 
0302 if sol.problem
0303     % Power flow equations may have a solution
0304     insolvable = 0;
0305 elseif sol.problem == 0
0306     % Power flow equations do not have a solution.
0307     insolvable = 1;
0308 end