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 %   $Id: insolvablepfsos_limitQ.m 2280 2014-01-17 23:28:37Z ray $
0039 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0040 %   and Ray Zimmerman, PSERC Cornell
0041 %   Copyright (c) 2013-2014 by Power System Engineering Research Center (PSERC)
0042 %
0043 %   This file is part of MATPOWER.
0044 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0045 %
0046 %   MATPOWER is free software: you can redistribute it and/or modify
0047 %   it under the terms of the GNU General Public License as published
0048 %   by the Free Software Foundation, either version 3 of the License,
0049 %   or (at your option) any later version.
0050 %
0051 %   MATPOWER is distributed in the hope that it will be useful,
0052 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0053 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0054 %   GNU General Public License for more details.
0055 %
0056 %   You should have received a copy of the GNU General Public License
0057 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0058 %
0059 %   Additional permission under GNU GPL version 3 section 7
0060 %
0061 %   If you modify MATPOWER, or any covered work, to interface with
0062 %   other modules (such as MATLAB code and MEX-files) available in a
0063 %   MATLAB(R) or comparable environment containing parts covered
0064 %   under other licensing terms, the licensors of MATPOWER grant
0065 %   you additional permission to convey the resulting work.
0066 
0067 define_constants;
0068 
0069 if nargin < 2
0070     mpopt = mpoption;
0071 end
0072 
0073 % Unpack options
0074 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0075 verbose             = mpopt.verbose;
0076 
0077 if ~have_fcn('yalmip')
0078     error('insolvablepfsos_limitQ: The software package YALMIP is required to run insolvablepfsos_limitQ. 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 %% Load data
0085 
0086 mpc = loadcase(mpc);
0087 mpc = ext2int(mpc);
0088 
0089 if toggle_dcline(mpc, 'status')
0090     error('insolvablepfsos_limitQ: DC lines are not implemented in insolvablepf');
0091 end
0092 
0093 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0094     warning('insolvablepfsos_limitQ: Angle difference constraints are not implemented. Ignoring angle difference constraints.');
0095 end
0096 
0097 nbus = size(mpc.bus,1);
0098 ngen = size(mpc.gen,1);
0099 
0100 % Some of the larger system (e.g., case2746wp) have generators
0101 % corresponding to buses that have bustype == PQ. Change these
0102 % to PV buses.
0103 for i=1:ngen
0104     busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0105     if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0106         mpc.bus(busidx,BUS_TYPE) = PV;
0107         if verbose >= 1
0108             warning('insolvablepfsos_limitQ: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0109         end
0110     end
0111 end
0112 
0113 % Buses may be listed as PV buses without associated generators. Change
0114 % these buses to PQ.
0115 for i=1:nbus
0116     if mpc.bus(i,BUS_TYPE) == PV
0117         genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0118         if isempty(genidx)
0119             mpc.bus(i,BUS_TYPE) = PQ;
0120             if verbose >= 1
0121                 warning('insolvablepfsos_limitQ: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0122             end
0123         end
0124     end
0125 end
0126 
0127 pq = find(mpc.bus(:,2) == PQ);
0128 pv = find(mpc.bus(:,2) == PV);
0129 ref = find(mpc.bus(:,2) == REF);
0130 
0131 refpv = sort([ref; pv]);
0132 pvpq = sort([pv; pq]);
0133 
0134 Y = makeYbus(mpc);
0135 G = real(Y);
0136 B = imag(Y);
0137 
0138 % Make vectors of power injections and voltage magnitudes
0139 S = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0140 P = real(S);
0141 Q = imag(S);
0142 
0143 Vmag = zeros(nbus,1);
0144 Vmag(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VG);
0145 
0146 % Determine the limits on the injected power at each bus
0147 Qmax = 1e3*ones(length(refpv),1);
0148 for i=1:length(refpv)
0149     genidx = find(mpc.gen(:,1) == refpv(i));
0150     Qmax(i) = (sum(mpc.gen(genidx,QMAX)) - mpc.bus(mpc.gen(genidx(1),1),QD)) / mpc.baseMVA;
0151 %     Qmin(i) = (sum(mpc.gen(genidx,QMIN)) - mpc.bus(mpc.gen(genidx(1),1),QD)) / mpc.baseMVA;
0152 end
0153 
0154 % bigM = 1e2;
0155 
0156 %% Create voltage variables
0157 % We need a Vd and Vq for each bus
0158 
0159 if verbose >= 2
0160     fprintf('Creating variables.\n');
0161 end
0162 
0163 Vd = sdpvar(nbus,1);
0164 Vq = sdpvar(nbus,1);
0165 
0166 % Variables for Q limits
0167 % Vplus2 = sdpvar(ngen,1);
0168 Vminus2 = sdpvar(ngen,1);
0169 x = sdpvar(ngen,1);
0170 
0171 % V = [1; Vd; Vq; Vplus2; Vminus2; x];
0172 V = [1; Vd; Vq; Vminus2; x];
0173 
0174 %% Create polynomials
0175 
0176 if verbose >= 2
0177     fprintf('Creating polynomials.\n');
0178 end
0179 
0180 deg = 0;
0181 sosdeg = 0;
0182 
0183 % Polynomials for active power
0184 for i=1:nbus-1
0185     [l,lc] = polynomial(V,deg,0);
0186     lambda(i) = l;
0187     clambda{i} = lc;
0188 end
0189 
0190 % Polynomials for reactive power
0191 for i=1:length(pq)
0192     [g,gc] = polynomial(V,deg,0);
0193     gamma(i) = g;
0194     cgamma{i} = gc;
0195 end
0196 
0197 % Polynomial for reference angle
0198 [delta, cdelta] = polynomial(V,deg,0);
0199 
0200 % Polynomials for Qlimits
0201 for i=1:length(refpv)
0202     [z,cz] = polynomial(V,deg,0);
0203     zeta{i}(1) = z;
0204     czeta{i}{1} = cz;
0205 
0206     [z,cz] = polynomial(V,deg,0);
0207     zeta{i}(2) = z;
0208     czeta{i}{2} = cz;
0209     
0210     [z,cz] = polynomial(V,deg,0);
0211     zeta{i}(3) = z;
0212     czeta{i}{3} = cz;
0213     
0214 %     [z,cz] = polynomial(V,deg,0);
0215 %     zeta{i}(4) = z;
0216 %     czeta{i}{4} = cz;
0217     
0218     [s1,c1] = polynomial(V,sosdeg,0);
0219     [s2,c2] = polynomial(V,sosdeg,0);
0220 %     [s3,c3] = polynomial(V,sosdeg,0);
0221 %     [s4,c4] = polynomial(V,sosdeg,0);
0222 %     [s5,c5] = polynomial(V,sosdeg,0);
0223     sigma{i}(1) = s1;
0224     csigma{i}{1} = c1;
0225     sigma{i}(2) = s2;
0226     csigma{i}{2} = c2;
0227 %     sigma{i}(3) = s3;
0228 %     csigma{i}{3} = c3;
0229 %     sigma{i}(4) = s4;
0230 %     csigma{i}{4} = c4;
0231 %     sigma{i}(5) = s5;
0232 %     csigma{i}{5} = c5;
0233 end
0234 
0235 %% Create power flow equation polynomials
0236 
0237 if verbose >= 2
0238     fprintf('Creating power flow equations.\n');
0239 end
0240 
0241 % Make P equations
0242 for k=1:nbus-1
0243     i = pvpq(k);
0244     
0245     % Take advantage of sparsity in forming polynomials
0246 %     fp(k) = Vd(i) * sum(G(:,i).*Vd - B(:,i).*Vq) + Vq(i) * sum(B(:,i).*Vd + G(:,i).*Vq);
0247     
0248     y = find(G(:,i) | B(:,i));
0249     for m=1:length(y)
0250         if exist('fp','var') && length(fp) == k
0251             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)));
0252         else
0253             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)));
0254         end
0255     end
0256 
0257 end
0258 
0259 % Make Q equations
0260 for k=1:length(pq)
0261     i = pq(k);
0262     
0263     % Take advantage of sparsity in forming polynomials
0264 %     fq(k) = Vd(i) * sum(-B(:,i).*Vd - G(:,i).*Vq) + Vq(i) * sum(G(:,i).*Vd - B(:,i).*Vq);
0265     
0266     y = find(G(:,i) | B(:,i));
0267     for m=1:length(y)
0268         if exist('fq','var') && length(fq) == k
0269             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)));
0270         else
0271             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)));
0272         end
0273     end
0274 end
0275 
0276 % Make V equations
0277 for k=1:length(refpv)
0278     i = refpv(k);
0279     fv(k) = Vd(i)^2 + Vq(i)^2;
0280 end
0281 
0282 % Generator reactive power limits
0283 for k=1:length(refpv)
0284     i = refpv(k);
0285     
0286     % Take advantage of sparsity in forming polynomials
0287 %     fq_gen(k) = Vd(i) * sum(-B(:,i).*Vd - G(:,i).*Vq) + Vq(i) * sum(G(:,i).*Vd - B(:,i).*Vq);
0288     
0289     y = find(G(:,i) | B(:,i));
0290     for m=1:length(y)
0291         if exist('fq_gen','var') && length(fq_gen) == k
0292             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)));
0293         else
0294             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)));
0295         end
0296     end
0297 end
0298 
0299 
0300 %% Form the test polynomial
0301 
0302 if verbose >= 2
0303     fprintf('Forming the test polynomial.\n');
0304 end
0305 
0306 sospoly = 1;
0307 
0308 % Add active power terms
0309 for i=1:nbus-1
0310     sospoly = sospoly + lambda(i) * (fp(i) - P(pvpq(i)));
0311 end
0312 
0313 % Add reactive power terms
0314 for i=1:length(pq)
0315     sospoly = sospoly + gamma(i) * (fq(i) - Q(pq(i)));
0316 end
0317 
0318 % Add angle reference term
0319 sospoly = sospoly + delta * Vq(ref);
0320 
0321 % Generator reactive power limits
0322 for k=1:length(refpv)
0323     i = refpv(k);
0324     
0325     % Both limits
0326 %     sospoly = sospoly + zeta{k}(1) * (Vmag(i)^2 + Vplus2(k) - Vminus2(k) - fv(k));
0327 %     sospoly = sospoly + zeta{k}(2) * (Vminus2(k) * x(k));
0328 %     sospoly = sospoly + zeta{k}(3) * (Vplus2(k) * (Qmax(k)-Qmin(k)-x(k)));
0329 %     sospoly = sospoly + zeta{k}(4) * (Qmax(k) - x(k) - fq_gen(k));
0330 %
0331 %     sospoly = sospoly + sigma{k}(1) * Vplus2(k);
0332 %     sospoly = sospoly + sigma{k}(2) * Vminus2(k);
0333 %     sospoly = sospoly + sigma{k}(3) * x(k);
0334 %     sospoly = sospoly + sigma{k}(4) * (Qmax(k) - Qmin(k) - x(k));
0335 %     sospoly = sospoly + sigma{k}(5) * (bigM - Vplus2(k));
0336     
0337     % Only upper reactive power limits
0338     sospoly = sospoly + zeta{k}(1) * (Vmag(i)^2 - Vminus2(k) - fv(k));
0339     sospoly = sospoly + zeta{k}(2) * (Vminus2(k) * x(k));
0340     sospoly = sospoly + zeta{k}(3) * (Qmax(k) - x(k) - fq_gen(k));
0341 
0342     sospoly = sospoly + sigma{k}(1) * Vminus2(k);
0343     sospoly = sospoly + sigma{k}(2) * x(k);
0344 end
0345 
0346 sospoly = -sospoly;
0347 
0348 %% Solve the SOS problem
0349 
0350 if verbose >= 2
0351     fprintf('Solving the sum of squares problem.\n');
0352 end
0353 
0354 sosopts = sdpsettings(sdpopts,'sos.newton',1,'sos.congruence',1);
0355 
0356 constraints = sos(sospoly);
0357 
0358 pvec = cdelta(:).';
0359 
0360 for i=1:length(lambda)
0361     pvec = [pvec clambda{i}(:).'];
0362 end
0363 
0364 for i=1:length(gamma)
0365     pvec = [pvec cgamma{i}(:).'];
0366 end
0367 
0368 for i=1:length(refpv)
0369     for k=1:length(czeta{i})
0370         pvec = [pvec czeta{i}{k}(:).'];
0371     end
0372     
0373     for k=1:length(csigma{i})
0374         pvec = [pvec csigma{i}{k}(:).'];
0375     end
0376     
0377     for k=1:length(csigma{i})
0378         constraints = [constraints; sos(sigma{i}(k))];
0379     end
0380 end
0381 
0382 % Preserve warning settings
0383 S = warning;
0384 
0385 sol = solvesos(constraints,[],sosopts,pvec);
0386 
0387 warning(S);
0388 
0389 if verbose >= 2
0390     fprintf('Solver exit message: %s\n',sol.info);
0391 end
0392 
0393 if sol.problem
0394     % Power flow equations may have a solution
0395     insolvable = 0;
0396 elseif sol.problem == 0
0397     % Power flow equations do not have a solution.
0398     insolvable = 1;
0399 end