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-2019, Power Systems Engineering Research Center (PSERC)
0039 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0040 %   and Ray Zimmerman, PSERC Cornell
0041 %
0042 %   This file is part of MATPOWER/mx-sdp_pf.
0043 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0044 %   See https://github.com/MATPOWER/mx-sdp_pf/ for more info.
0045 
0046 define_constants;
0047 
0048 if nargin < 2
0049     mpopt = mpoption;
0050 end
0051 
0052 % Unpack options
0053 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0054 verbose             = mpopt.verbose;
0055 
0056 if ~have_fcn('yalmip')
0057     error('insolvablepfsos_limitQ: The software package YALMIP is required to run insolvablepfsos_limitQ. 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 %% Load data
0064 
0065 mpc = loadcase(mpc);
0066 mpc = ext2int(mpc);
0067 
0068 if toggle_dcline(mpc, 'status')
0069     error('insolvablepfsos_limitQ: DC lines are not implemented in insolvablepf');
0070 end
0071 
0072 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0073     warning('insolvablepfsos_limitQ: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0074 end
0075 
0076 nbus = size(mpc.bus,1);
0077 ngen = size(mpc.gen,1);
0078 
0079 % Some of the larger system (e.g., case2746wp) have generators
0080 % corresponding to buses that have bustype == PQ. Change these
0081 % to PV buses.
0082 for i=1:ngen
0083     busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0084     if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0085         mpc.bus(busidx,BUS_TYPE) = PV;
0086         if verbose >= 1
0087             warning('insolvablepfsos_limitQ: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0088         end
0089     end
0090 end
0091 
0092 % Buses may be listed as PV buses without associated generators. Change
0093 % these buses to PQ.
0094 for i=1:nbus
0095     if mpc.bus(i,BUS_TYPE) == PV
0096         genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0097         if isempty(genidx)
0098             mpc.bus(i,BUS_TYPE) = PQ;
0099             if verbose >= 1
0100                 warning('insolvablepfsos_limitQ: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0101             end
0102         end
0103     end
0104 end
0105 
0106 pq = find(mpc.bus(:,2) == PQ);
0107 pv = find(mpc.bus(:,2) == PV);
0108 ref = find(mpc.bus(:,2) == REF);
0109 
0110 refpv = sort([ref; pv]);
0111 pvpq = sort([pv; pq]);
0112 
0113 Y = makeYbus(mpc);
0114 G = real(Y);
0115 B = imag(Y);
0116 
0117 % Make vectors of power injections and voltage magnitudes
0118 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0119 P = real(S);
0120 Q = imag(S);
0121 
0122 Vmag = zeros(nbus,1);
0123 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0124 
0125 % Determine the limits on the injected power at each bus
0126 Qmax = 1e3*ones(length(refpv),1);
0127 for i=1:length(refpv)
0128     genidx = find(mpc.gen(:,1) == refpv(i));
0129     Qmax(i) = (sum(mpc.gen(genidx,QMAX)) - mpc.bus(mpc.gen(genidx(1),1),QD)) / mpc.baseMVA;
0130 %     Qmin(i) = (sum(mpc.gen(genidx,QMIN)) - mpc.bus(mpc.gen(genidx(1),1),QD)) / mpc.baseMVA;
0131 end
0132 
0133 % bigM = 1e2;
0134 
0135 %% Create voltage variables
0136 % We need a Vd and Vq for each bus
0137 
0138 if verbose >= 2
0139     fprintf('Creating variables.\n');
0140 end
0141 
0142 Vd = sdpvar(nbus,1);
0143 Vq = sdpvar(nbus,1);
0144 
0145 % Variables for Q limits
0146 % Vplus2 = sdpvar(ngen,1);
0147 Vminus2 = sdpvar(ngen,1);
0148 x = sdpvar(ngen,1);
0149 
0150 % V = [1; Vd; Vq; Vplus2; Vminus2; x];
0151 V = [1; Vd; Vq; Vminus2; x];
0152 
0153 %% Create polynomials
0154 
0155 if verbose >= 2
0156     fprintf('Creating polynomials.\n');
0157 end
0158 
0159 deg = 0;
0160 sosdeg = 0;
0161 
0162 % Polynomials for active power
0163 for i=1:nbus-1
0164     [l,lc] = polynomial(V,deg,0);
0165     lambda(i) = l;
0166     clambda{i} = lc;
0167 end
0168 
0169 % Polynomials for reactive power
0170 for i=1:length(pq)
0171     [g,gc] = polynomial(V,deg,0);
0172     gamma(i) = g;
0173     cgamma{i} = gc;
0174 end
0175 
0176 % Polynomial for reference angle
0177 [delta, cdelta] = polynomial(V,deg,0);
0178 
0179 % Polynomials for Qlimits
0180 for i=1:length(refpv)
0181     [z,cz] = polynomial(V,deg,0);
0182     zeta{i} = z;    %% RDZ: required by Octave (4.0.3) to initialize properly
0183     zeta{i}(1) = z;
0184     czeta{i}{1} = cz;
0185 
0186     [z,cz] = polynomial(V,deg,0);
0187     zeta{i}(2) = z;
0188     czeta{i}{2} = cz;
0189     
0190     [z,cz] = polynomial(V,deg,0);
0191     zeta{i}(3) = z;
0192     czeta{i}{3} = cz;
0193     
0194 %     [z,cz] = polynomial(V,deg,0);
0195 %     zeta{i}(4) = z;
0196 %     czeta{i}{4} = cz;
0197     
0198     [s1,c1] = polynomial(V,sosdeg,0);
0199     [s2,c2] = polynomial(V,sosdeg,0);
0200 %     [s3,c3] = polynomial(V,sosdeg,0);
0201 %     [s4,c4] = polynomial(V,sosdeg,0);
0202 %     [s5,c5] = polynomial(V,sosdeg,0);
0203     sigma{i} = s1;  %% RDZ: required by Octave (4.0.3) to initialize properly
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 if have_fcn('matlab', 'vnum') >= 8.006 && have_fcn('cplex') && ...
0366         have_fcn('cplex', 'vnum') <= 12.006003
0367     warning('OFF', 'MATLAB:lang:badlyScopedReturnValue');
0368 end
0369 
0370 sol = solvesos(constraints,[],sosopts,pvec);
0371 
0372 if ~have_fcn('octave') || have_fcn('octave', 'vnum') >= 4.001
0373     %% (avoid bug in Octave 4.0.x, where warning state is left corrupted)
0374     warning(S);
0375 end
0376 
0377 if verbose >= 2
0378     fprintf('Solver exit message: %s\n',sol.info);
0379 end
0380 
0381 if sol.problem
0382     % Power flow equations may have a solution
0383     insolvable = 0;
0384 elseif sol.problem == 0
0385     % Power flow equations do not have a solution.
0386     insolvable = 1;
0387 end