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 %   $Id: insolvablepfsos.m 2280 2014-01-17 23:28:37Z ray $
0035 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0036 %   Copyright (c) 2013-2014 by Power System Engineering Research Center (PSERC)
0037 %
0038 %   This file is part of MATPOWER.
0039 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0040 %
0041 %   MATPOWER is free software: you can redistribute it and/or modify
0042 %   it under the terms of the GNU General Public License as published
0043 %   by the Free Software Foundation, either version 3 of the License,
0044 %   or (at your option) any later version.
0045 %
0046 %   MATPOWER is distributed in the hope that it will be useful,
0047 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0048 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0049 %   GNU General Public License for more details.
0050 %
0051 %   You should have received a copy of the GNU General Public License
0052 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0053 %
0054 %   Additional permission under GNU GPL version 3 section 7
0055 %
0056 %   If you modify MATPOWER, or any covered work, to interface with
0057 %   other modules (such as MATLAB code and MEX-files) available in a
0058 %   MATLAB(R) or comparable environment containing parts covered
0059 %   under other licensing terms, the licensors of MATPOWER grant
0060 %   you additional permission to convey the resulting work.
0061 
0062 define_constants;
0063 
0064 if nargin < 2
0065     mpopt = mpoption;
0066 end
0067 
0068 % Unpack options
0069 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0070 verbose             = mpopt.verbose;
0071 enforce_Qlimits     = mpopt.pf.enforce_q_lims;
0072 
0073 if enforce_Qlimits > 0
0074     enforce_Qlimits = 1;
0075 end
0076 
0077 if ~have_fcn('yalmip')
0078     error('insolvablepfsos: The software package YALMIP is required to run insolvablepfsos. See http://users.isy.liu.se/johanl/yalmip/');
0079 end
0080 
0081 % set YALMIP options struct in SDP_PF (for details, see help sdpsettings)
0082 sdpopts = yalmip_options([], mpopt);
0083 
0084 %% Handle generator reactive power limits
0085 % If no generator reactive power limits are specified, use this code
0086 % directly. If generator reactive power limits are to be enforced, use a
0087 % function that considers reactive power limited generators.
0088 
0089 if enforce_Qlimits
0090     if verbose > 0
0091         fprintf('Generator reactive power limits are enforced. Using function insolvablepfsos_limitQ.\n');
0092     end
0093     insolvable = insolvablepfsos_limitQ(mpc,mpopt);
0094     return;
0095 end
0096 
0097 %% Load data
0098 
0099 mpc = loadcase(mpc);
0100 mpc = ext2int(mpc);
0101 
0102 if toggle_dcline(mpc, 'status')
0103     error('insolvablepfsos: DC lines are not implemented in insolvablepf');
0104 end
0105 
0106 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0107     warning('insolvablepfsos: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0108 end
0109 
0110 nbus = size(mpc.bus,1);
0111 ngen = size(mpc.gen,1);
0112 
0113 % Some of the larger system (e.g., case2746wp) have generators
0114 % corresponding to buses that have bustype == PQ. Change these
0115 % to PV buses.
0116 for i=1:ngen
0117     busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0118     if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0119         mpc.bus(busidx,BUS_TYPE) = PV;
0120         if verbose >= 1
0121             warning('insolvablepfsos: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0122         end
0123     end
0124 end
0125 
0126 % Buses may be listed as PV buses without associated generators. Change
0127 % these buses to PQ.
0128 for i=1:nbus
0129     if mpc.bus(i,BUS_TYPE) == PV
0130         genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0131         if isempty(genidx)
0132             mpc.bus(i,BUS_TYPE) = PQ;
0133             if verbose >= 1
0134                 warning('insolvablepfsos: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0135             end
0136         end
0137     end
0138 end
0139 
0140 pq = find(mpc.bus(:,2) == PQ);
0141 pv = find(mpc.bus(:,2) == PV);
0142 ref = find(mpc.bus(:,2) == REF);
0143 
0144 refpv = sort([ref; pv]);
0145 pvpq = sort([pv; pq]);
0146 
0147 Y = makeYbus(mpc);
0148 G = real(Y);
0149 B = imag(Y);
0150 
0151 % Make vectors of power injections and voltage magnitudes
0152 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0153 P = real(S);
0154 Q = imag(S);
0155 
0156 Vmag = zeros(nbus,1);
0157 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0158 
0159 
0160 %% Create voltage variables
0161 % We need a Vd and Vq for each bus
0162 
0163 if verbose >= 2
0164     fprintf('Creating variables.\n');
0165 end
0166 
0167 Vd = sdpvar(nbus,1);
0168 Vq = sdpvar(nbus,1);
0169 
0170 V = [1; Vd; Vq];
0171 
0172 %% Create polynomials
0173 
0174 if verbose >= 2
0175     fprintf('Creating polynomials.\n');
0176 end
0177 
0178 deg = 0;
0179 
0180 % Polynomials for active power
0181 for i=1:nbus-1
0182     [l,lc] = polynomial(V,deg,0);
0183     lambda(i) = l;
0184     clambda{i} = lc;
0185 end
0186 
0187 % Polynomials for reactive power
0188 for i=1:length(pq)
0189     [g,gc] = polynomial(V,deg,0);
0190     gamma(i) = g;
0191     cgamma{i} = gc;
0192 end
0193 
0194 % Polynomials for voltage magnitude
0195 for i=1:length(refpv)
0196     [m,mc] = polynomial(V,deg,0);
0197     mu(i) = m;
0198     cmu{i} = mc;
0199 end
0200 
0201 % Polynomial for reference angle
0202 [delta, cdelta] = polynomial(V,deg,0);
0203 
0204 
0205 %% Create power flow equation polynomials
0206 
0207 if verbose >= 2
0208     fprintf('Creating power flow equations.\n');
0209 end
0210 
0211 % Make P equations
0212 for k=1:nbus-1
0213     i = pvpq(k);
0214     
0215     % Take advantage of sparsity in forming polynomials
0216 %     fp(k) = Vd(i) * sum(G(:,i).*Vd - B(:,i).*Vq) + Vq(i) * sum(B(:,i).*Vd + G(:,i).*Vq);
0217     
0218     y = find(G(:,i) | B(:,i));
0219     for m=1:length(y)
0220         if exist('fp','var') && length(fp) == k
0221             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)));
0222         else
0223             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)));
0224         end
0225     end
0226 
0227 end
0228 
0229 % Make Q equations
0230 for k=1:length(pq)
0231     i = pq(k);
0232     
0233     % Take advantage of sparsity in forming polynomials
0234 %     fq(k) = Vd(i) * sum(-B(:,i).*Vd - G(:,i).*Vq) + Vq(i) * sum(G(:,i).*Vd - B(:,i).*Vq);
0235     
0236     y = find(G(:,i) | B(:,i));
0237     for m=1:length(y)
0238         if exist('fq','var') && length(fq) == k
0239             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)));
0240         else
0241             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)));
0242         end
0243     end
0244 end
0245 
0246 % Make V equations
0247 for k=1:length(refpv)
0248     i = refpv(k);
0249     fv(k) = Vd(i)^2 + Vq(i)^2;
0250 end
0251 
0252 
0253 %% Form the test polynomial
0254 
0255 if verbose >= 2
0256     fprintf('Forming the test polynomial.\n');
0257 end
0258 
0259 sospoly = 1;
0260 
0261 % Add active power terms
0262 for i=1:nbus-1
0263     sospoly = sospoly + lambda(i) * (fp(i) - P(pvpq(i)));
0264 end
0265 
0266 % Add reactive power terms
0267 for i=1:length(pq)
0268     sospoly = sospoly + gamma(i) * (fq(i) - Q(pq(i)));
0269 end
0270 
0271 % Add voltage magnitude terms
0272 for i=1:length(refpv)
0273     sospoly = sospoly + mu(i) * (fv(i) - Vmag(refpv(i))^2);
0274 end
0275 
0276 % Add angle reference term
0277 sospoly = sospoly + delta * Vq(ref);
0278 
0279 sospoly = -sospoly;
0280 
0281 %% Solve the SOS problem
0282 
0283 if verbose >= 2
0284     fprintf('Solving the sum of squares problem.\n');
0285 end
0286 
0287 sosopts = sdpsettings(sdpopts,'sos.newton',1,'sos.congruence',1);
0288 
0289 constraints = sos(sospoly);
0290 
0291 pvec = cdelta(:).';
0292 
0293 for i=1:length(lambda)
0294     pvec = [pvec clambda{i}(:).'];
0295 end
0296 
0297 for i=1:length(gamma)
0298     pvec = [pvec cgamma{i}(:).'];
0299 end
0300 
0301 for i=1:length(mu)
0302     pvec = [pvec cmu{i}(:).'];
0303 end
0304 
0305 % Preserve warning settings
0306 S = warning;
0307 
0308 sol = solvesos(constraints,[],sosopts,pvec);
0309 
0310 warning(S);
0311 
0312 if verbose >= 2
0313     fprintf('Solver exit message: %s\n',sol.info);
0314 end
0315 
0316 if sol.problem
0317     % Power flow equations may have a solution
0318     insolvable = 0;
0319 elseif sol.problem == 0
0320     % Power flow equations do not have a solution.
0321     insolvable = 1;
0322 end