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-2016, 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.
0043 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0044 %   See http://www.pserc.cornell.edu/matpower/ 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 http://users.isy.liu.se/johanl/yalmip/');
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}(1) = z;
0183     czeta{i}{1} = cz;
0184 
0185     [z,cz] = polynomial(V,deg,0);
0186     zeta{i}(2) = z;
0187     czeta{i}{2} = cz;
0188     
0189     [z,cz] = polynomial(V,deg,0);
0190     zeta{i}(3) = z;
0191     czeta{i}{3} = cz;
0192     
0193 %     [z,cz] = polynomial(V,deg,0);
0194 %     zeta{i}(4) = z;
0195 %     czeta{i}{4} = cz;
0196     
0197     [s1,c1] = polynomial(V,sosdeg,0);
0198     [s2,c2] = polynomial(V,sosdeg,0);
0199 %     [s3,c3] = polynomial(V,sosdeg,0);
0200 %     [s4,c4] = polynomial(V,sosdeg,0);
0201 %     [s5,c5] = polynomial(V,sosdeg,0);
0202     sigma{i}(1) = s1;
0203     csigma{i}{1} = c1;
0204     sigma{i}(2) = s2;
0205     csigma{i}{2} = c2;
0206 %     sigma{i}(3) = s3;
0207 %     csigma{i}{3} = c3;
0208 %     sigma{i}(4) = s4;
0209 %     csigma{i}{4} = c4;
0210 %     sigma{i}(5) = s5;
0211 %     csigma{i}{5} = c5;
0212 end
0213 
0214 %% Create power flow equation polynomials
0215 
0216 if verbose >= 2
0217     fprintf('Creating power flow equations.\n');
0218 end
0219 
0220 % Make P equations
0221 for k=1:nbus-1
0222     i = pvpq(k);
0223     
0224     % Take advantage of sparsity in forming polynomials
0225 %     fp(k) = Vd(i) * sum(G(:,i).*Vd - B(:,i).*Vq) + Vq(i) * sum(B(:,i).*Vd + G(:,i).*Vq);
0226     
0227     y = find(G(:,i) | B(:,i));
0228     for m=1:length(y)
0229         if exist('fp','var') && length(fp) == k
0230             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)));
0231         else
0232             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         end
0234     end
0235 
0236 end
0237 
0238 % Make Q equations
0239 for k=1:length(pq)
0240     i = pq(k);
0241     
0242     % Take advantage of sparsity in forming polynomials
0243 %     fq(k) = Vd(i) * sum(-B(:,i).*Vd - G(:,i).*Vq) + Vq(i) * sum(G(:,i).*Vd - B(:,i).*Vq);
0244     
0245     y = find(G(:,i) | B(:,i));
0246     for m=1:length(y)
0247         if exist('fq','var') && length(fq) == k
0248             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)));
0249         else
0250             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         end
0252     end
0253 end
0254 
0255 % Make V equations
0256 for k=1:length(refpv)
0257     i = refpv(k);
0258     fv(k) = Vd(i)^2 + Vq(i)^2;
0259 end
0260 
0261 % Generator reactive power limits
0262 for k=1:length(refpv)
0263     i = refpv(k);
0264     
0265     % Take advantage of sparsity in forming polynomials
0266 %     fq_gen(k) = Vd(i) * sum(-B(:,i).*Vd - G(:,i).*Vq) + Vq(i) * sum(G(:,i).*Vd - B(:,i).*Vq);
0267     
0268     y = find(G(:,i) | B(:,i));
0269     for m=1:length(y)
0270         if exist('fq_gen','var') && length(fq_gen) == k
0271             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)));
0272         else
0273             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         end
0275     end
0276 end
0277 
0278 
0279 %% Form the test polynomial
0280 
0281 if verbose >= 2
0282     fprintf('Forming the test polynomial.\n');
0283 end
0284 
0285 sospoly = 1;
0286 
0287 % Add active power terms
0288 for i=1:nbus-1
0289     sospoly = sospoly + lambda(i) * (fp(i) - P(pvpq(i)));
0290 end
0291 
0292 % Add reactive power terms
0293 for i=1:length(pq)
0294     sospoly = sospoly + gamma(i) * (fq(i) - Q(pq(i)));
0295 end
0296 
0297 % Add angle reference term
0298 sospoly = sospoly + delta * Vq(ref);
0299 
0300 % Generator reactive power limits
0301 for k=1:length(refpv)
0302     i = refpv(k);
0303     
0304     % Both limits
0305 %     sospoly = sospoly + zeta{k}(1) * (Vmag(i)^2 + Vplus2(k) - Vminus2(k) - fv(k));
0306 %     sospoly = sospoly + zeta{k}(2) * (Vminus2(k) * x(k));
0307 %     sospoly = sospoly + zeta{k}(3) * (Vplus2(k) * (Qmax(k)-Qmin(k)-x(k)));
0308 %     sospoly = sospoly + zeta{k}(4) * (Qmax(k) - x(k) - fq_gen(k));
0309 %
0310 %     sospoly = sospoly + sigma{k}(1) * Vplus2(k);
0311 %     sospoly = sospoly + sigma{k}(2) * Vminus2(k);
0312 %     sospoly = sospoly + sigma{k}(3) * x(k);
0313 %     sospoly = sospoly + sigma{k}(4) * (Qmax(k) - Qmin(k) - x(k));
0314 %     sospoly = sospoly + sigma{k}(5) * (bigM - Vplus2(k));
0315     
0316     % Only upper reactive power limits
0317     sospoly = sospoly + zeta{k}(1) * (Vmag(i)^2 - Vminus2(k) - fv(k));
0318     sospoly = sospoly + zeta{k}(2) * (Vminus2(k) * x(k));
0319     sospoly = sospoly + zeta{k}(3) * (Qmax(k) - x(k) - fq_gen(k));
0320 
0321     sospoly = sospoly + sigma{k}(1) * Vminus2(k);
0322     sospoly = sospoly + sigma{k}(2) * x(k);
0323 end
0324 
0325 sospoly = -sospoly;
0326 
0327 %% Solve the SOS problem
0328 
0329 if verbose >= 2
0330     fprintf('Solving the sum of squares problem.\n');
0331 end
0332 
0333 sosopts = sdpsettings(sdpopts,'sos.newton',1,'sos.congruence',1);
0334 
0335 constraints = sos(sospoly);
0336 
0337 pvec = cdelta(:).';
0338 
0339 for i=1:length(lambda)
0340     pvec = [pvec clambda{i}(:).'];
0341 end
0342 
0343 for i=1:length(gamma)
0344     pvec = [pvec cgamma{i}(:).'];
0345 end
0346 
0347 for i=1:length(refpv)
0348     for k=1:length(czeta{i})
0349         pvec = [pvec czeta{i}{k}(:).'];
0350     end
0351     
0352     for k=1:length(csigma{i})
0353         pvec = [pvec csigma{i}{k}(:).'];
0354     end
0355     
0356     for k=1:length(csigma{i})
0357         constraints = [constraints; sos(sigma{i}(k))];
0358     end
0359 end
0360 
0361 % Preserve warning settings
0362 S = warning;
0363 if have_fcn('matlab', 'vnum') >= 8.006 && have_fcn('cplex') && ...
0364         have_fcn('cplex', 'vnum') <= 12.006003
0365     warning('OFF', 'MATLAB:lang:badlyScopedReturnValue');
0366 end
0367 
0368 sol = solvesos(constraints,[],sosopts,pvec);
0369 
0370 warning(S);
0371 
0372 if verbose >= 2
0373     fprintf('Solver exit message: %s\n',sol.info);
0374 end
0375 
0376 if sol.problem
0377     % Power flow equations may have a solution
0378     insolvable = 0;
0379 elseif sol.problem == 0
0380     % Power flow equations do not have a solution.
0381     insolvable = 1;
0382 end