insolvablepfsos_limitQ

INSOLVABLEPFSOS_LIMITQ A sufficient condition for power flow insolvability
considering generator reactive power limits using sum of squares programming

   [INSOLVABLE] = INSOLVABLEPFSOS_LIMITQ(MPC,MPOPT)

   Uses sum of squares programming to generate an infeasibility 
   certificate for the power flow equations considering reactive power
   limited generators. 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.

   Note that only upper limits on reactive power generation are currently
   considered.

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