Description of insolvablepf

0001 function [insolvable,Vslack_min,sigma,eta,mineigratio] = insolvablepf(mpc,mpopt)
0002 %INSOLVABLEPF A sufficient condition for power flow insolvability
0003 %
0004 %   [INSOLVABLE,VSLACK_MIN,SIGMA,ETA,MINEIGRATIO] = INSOLVABLEPF(MPC,MPOPT)
0005 %
0006 %   Evaluates a sufficient condition for insolvability of the power flow
0007 %   equations. The voltage at the slack bus is minimized (with
0008 %   proportional changes in voltage magnitudes at PV buses) using a
0009 %   semidefinite programming relaxation of the power flow equations. If the
0010 %   minimum achievable slack bus voltage is greater than the specified slack
0011 %   bus voltage, no power flow solutions exist. The converse does not
0012 %   necessarily hold; a power flow solution may not exist for cases
0013 %   where the output insolvable is equal to 0. See [1] and [2] for further
0014 %   details.
0015 %
0016 %   Inputs:
0017 %       MPC : A MATPOWER case specifying the desired power flow equations.
0018 %       MPOPT : A MATPOWER options struct. If not specified, it is
0019 %           assumed to be the default mpoption.
0020 %
0021 %   Outputs:
0022 %       INSOLVABLE : Binary variable. A value of 1 indicates that the
0023 %           specified power flow equations are insolvable, while a value of
0024 %           0 means that the insolvability condition is indeterminant (a
0025 %           solution may or may not exist).
0026 %       VSLACK_MIN : Minimum possible slack voltage obtained from the
0027 %           semidefinite programming relaxation. The power flow equations
0028 %           are insolvable if Vslack_min > V0, where V0 is the specified
0029 %           voltage at the slack bus.
0030 %       SIGMA : Controlled voltage margin to the power flow solvability
0031 %           boundary.
0032 %       ETA : Power injection margin to the power flow solvability
0033 %           boundary in the profile of a uniform, constant power factor
0034 %           change in power injections.
0035 %       MINEIGRATIO : A measure of satisfaction of the rank relaxation.
0036 %           Large values indicate satisfaction. (Note that satisfaction of
0037 %           the rank relaxation is not required for correctness of the
0038 %           insolvability condition).
0039 %
0040 %   Note that this function uses a matrix completion decomposition and is
0041 %   therefore suitable for large systems.
0042 %
0043 % [1] D.K. Molzahn, B.C. Lesieutre, and C.L. DeMarco, "A Sufficient Condition
0044 %     for Power Flow Insolvability with Applications to Voltage Stability
0045 %     Margins," IEEE Transactions on Power Systems, vol. 28, no. 3,
0046 %     pp. 2592-2601, August 2013.
0047 %
0048 % [2] D.K. Molzahn, B.C. Lesieutre, and C.L. DeMarco, "A Sufficient Condition
0049 %     for Power Flow Insolvability with Applications to Voltage Stability
0050 %     Margins," University of Wisconsin-Madison Department of Electrical
0051 %     and Computer Engineering, Tech. Rep. ECE-12-01, 2012, [Online].
0052 %     Available: https://arxiv.org/abs/1204.6285.
0053 
0054 %   MATPOWER
0055 %   Copyright (c) 2013-2019, Power Systems Engineering Research Center (PSERC)
0056 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0057 %   and Ray Zimmerman, PSERC Cornell
0058 %
0059 %   This file is part of MATPOWER/mx-sdp_pf.
0060 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0061 %   See https://github.com/MATPOWER/mx-sdp_pf/ for more info.
0062 
0063 if nargin < 2
0064     mpopt = mpoption;
0065 end
0066 
0067 %% define undocumented MATLAB function ismembc() if not available (e.g. Octave)
0068 if exist('ismembc')
0069     ismembc_ = @ismembc;
0070 else
0071     ismembc_ = @ismembc_octave;
0072 end
0073 
0074 mpc = loadcase(mpc);
0075 mpc = ext2int(mpc);
0076 
0077 % Unpack options
0078 ignore_angle_lim    = mpopt.opf.ignore_angle_lim;
0079 verbose             = mpopt.verbose;
0080 enforce_Qlimits     = mpopt.pf.enforce_q_lims;
0081 maxNumberOfCliques  = mpopt.sdp_pf.max_number_of_cliques;   %% Maximum number of maximal cliques
0082 ndisplay            = mpopt.sdp_pf.ndisplay;        %% Determine display frequency of diagonastic information
0083 cholopts.dense      = mpopt.sdp_pf.choldense;       %% Cholesky factorization options
0084 cholopts.aggressive = mpopt.sdp_pf.cholaggressive;  %% Cholesky factorization options
0085 
0086 if enforce_Qlimits > 0
0087     enforce_Qlimits = 1;
0088 end
0089 
0090 if ~have_feature('yalmip')
0091     error('insolvablepf: The software package YALMIP is required to run insolvablepf. See https://yalmip.github.io');
0092 end
0093 
0094 % set YALMIP options struct in SDP_PF (for details, see help sdpsettings)
0095 sdpopts = yalmip_options([], mpopt);
0096 
0097 %% Handle generator reactive power limits
0098 % If no generator reactive power limits are specified, use this code
0099 % directly. If generator reactive power limits are to be enforced, use the
0100 % mixed integer semidefinite programming code. This code is only applicable
0101 % for small systems (the 57-bus system is really pushing the limits)
0102 
0103 if enforce_Qlimits
0104     if verbose > 0
0105         fprintf('Generator reactive power limits are enforced. Using function insolvablepf_limitQ.\n');
0106     end
0107     [insolvable,eta,mineigratio] = insolvablepf_limitQ(mpc,mpopt);
0108     Vslack_min = nan;
0109     sigma = nan;
0110     return;
0111 end
0112 
0113 
0114 %% define named indices into data matrices
0115 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0116     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0117 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0118     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0119     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0120 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0121     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0122     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0123 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0124 
0125 
0126 %% Load mpc data
0127 if isfield(mpc, 'userfcn') && length(mpc.userfcn) > 0 && ...
0128         isfield(mpc.userfcn, 'formulation')
0129     c = cellfun(@func2str, {mpc.userfcn.formulation.fcn}, 'UniformOutput', 0);
0130     if strfind(strcat(c{:}), 'userfcn_dcline_formulation')
0131         error('insolvablepf: DC lines are not implemented in insolvablepf');
0132     end
0133 end
0134 
0135 if toggle_dcline(mpc, 'status')
0136     error('insolvablepf: DC lines are not implemented in insolvablepf');
0137 end
0138 
0139 nbus = size(mpc.bus,1);
0140 ngen = size(mpc.gen,1);
0141 nbranch = size(mpc.branch,1);
0142 
0143 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0144     warning('insolvablepf: Angle difference constraints are not implemented in SDP_PF. Ignoring angle difference constraints.');
0145 end
0146 
0147 % Some of the larger system (e.g., case2746wp) have generators
0148 % corresponding to buses that have bustype == PQ. Change these
0149 % to PV buses.
0150 for i=1:ngen
0151     busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0152     if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0153         mpc.bus(busidx,BUS_TYPE) = PV;
0154         if verbose >= 1
0155             warning('insolvablepf: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0156         end
0157     end
0158 end
0159 
0160 % Buses may be listed as PV buses without associated generators. Change
0161 % these buses to PQ.
0162 for i=1:nbus
0163     if mpc.bus(i,BUS_TYPE) == PV
0164         genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0165         if isempty(genidx)
0166             mpc.bus(i,BUS_TYPE) = PQ;
0167             if verbose >= 1
0168                 warning('insolvablepf: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0169             end
0170         end
0171     end
0172 end
0173 
0174 
0175 
0176 %% Determine Maximal Cliques
0177 %% Step 1: Cholesky factorization to obtain chordal extension
0178 % Use a minimum degree permutation to obtain a sparse chordal extension.
0179 
0180 if maxNumberOfCliques ~= 1
0181     [Ainc] = makeIncidence(mpc.bus,mpc.branch);
0182     sparsity_pattern = abs(Ainc.'*Ainc) + eye(nbus,nbus);
0183     per = amd(sparsity_pattern,cholopts);
0184     [Rchol,p] = chol(sparsity_pattern(per,per),'lower');
0185 
0186     if p ~= 0
0187         error('insolvablepf: sparsity_pattern not positive definite!');
0188     end
0189 else
0190     per = 1:nbus;
0191 end
0192 
0193 % Rearrange mpc to the same order as the permutation per
0194 mpc.bus = mpc.bus(per,:);
0195 mpc.bus(:,BUS_I) = 1:nbus;
0196 
0197 for i=1:ngen
0198     mpc.gen(i,GEN_BUS) = find(per == mpc.gen(i,GEN_BUS));
0199 end
0200 
0201 [mpc.gen genidx] = sortrows(mpc.gen,1);
0202 mpc.gencost = mpc.gencost(genidx,:);
0203 
0204 for i=1:nbranch
0205     mpc.branch(i,F_BUS) = find(per == mpc.branch(i,F_BUS));
0206     mpc.branch(i,T_BUS) = find(per == mpc.branch(i,T_BUS));
0207 end
0208 % -------------------------------------------------------------------------
0209 
0210 Sbase = mpc.baseMVA;
0211 
0212 % Create vectors of power injections and voltage magnitudes
0213 Qinj = -mpc.bus(:,QD) / Sbase;
0214 Vmag = mpc.bus(:,VM);
0215 
0216 Pd = mpc.bus(:,PD) / Sbase;
0217 Pg = zeros(nbus,1);
0218 for i=1:nbus
0219     genidx = find(mpc.gen(:,GEN_BUS) == i);
0220     if ~isempty(genidx)
0221         Pg(i) = sum(mpc.gen(genidx,PG)) / Sbase;
0222         Vmag(i) = mpc.gen(genidx(1),VG);
0223     end
0224 end
0225 Pinj = Pg - Pd;
0226 
0227 slackbus_idx = find(mpc.bus(:,BUS_TYPE) == 3);
0228 Vslack = Vmag(slackbus_idx);
0229 
0230 %% Step 2: Compute maximal cliques of chordal extension
0231 % Use adjacency matrix made from Rchol
0232 
0233 if maxNumberOfCliques ~= 1 && nbus > 3
0234     
0235     % Build adjacency matrix
0236     [f,t] = find(Rchol - diag(diag(Rchol)));
0237     Aadj = sparse(f,t,ones(length(f),1),nbus,nbus);
0238     Aadj = max(Aadj,Aadj');
0239 
0240     % Determine maximal cliques
0241     [MC,ischordal] = maxCardSearch(Aadj);
0242     
0243     if ~ischordal
0244         % This should never happen since the Cholesky decomposition should
0245         % always yield a chordal graph.
0246         error('insolvablepf: Chordal extension adjacency matrix is not chordal!');
0247     end
0248     
0249     for i=1:size(MC,2)
0250        maxclique{i} = find(MC(:,i));
0251     end
0252 
0253 else
0254     maxclique{1} = 1:nbus;
0255     nmaxclique = 1;
0256     E = [];
0257 end
0258 
0259 
0260 %% Step 3: Make clique tree and combine maximal cliques
0261 
0262 if maxNumberOfCliques ~= 1 && nbus > 3
0263     % Create a graph of maximal cliques, where cost between each maximal clique
0264     % is the number of shared buses in the corresponding maximal cliques.
0265     nmaxclique = length(maxclique);
0266     cliqueCost = sparse(nmaxclique,nmaxclique);
0267     for i=1:nmaxclique
0268         maxcliquei = maxclique{i};
0269         for k=i+1:nmaxclique
0270 
0271             cliqueCost(i,k) = sum(ismembc_(maxcliquei,maxclique{k}));
0272 
0273             % Slower alternative that doesn't use undocumented MATLAB function
0274             % cliqueCost(i,k) = length(intersect(maxcliquei,maxclique{k}));
0275 
0276         end
0277     end
0278     cliqueCost = max(cliqueCost,cliqueCost.');
0279 
0280     % Calculate the maximal spanning tree
0281     cliqueCost = -cliqueCost;
0282     [E] = prim(cliqueCost);
0283 
0284     % Combine maximal cliques
0285     if verbose >= 2
0286         [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,ndisplay);
0287     else
0288         [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,inf);
0289     end
0290     nmaxclique = length(maxclique);
0291 end
0292 
0293 %% Create SDP relaxation
0294 
0295 tic;
0296 
0297 constraints = [];
0298 
0299 % Control growth of Wref variables
0300 dd_blk_size = 50*nbus;
0301 dq_blk_size = 2*dd_blk_size;
0302 
0303 Wref_dd = zeros(dd_blk_size,2); % For terms like Vd1*Vd1 and Vd1*Vd2
0304 lWref_dd = dd_blk_size;
0305 matidx_dd = zeros(dd_blk_size,3);
0306 dd_ptr = 0;
0307 
0308 Wref_dq = zeros(dq_blk_size,2); % For terms like Vd1*Vd1 and Vd1*Vd2
0309 lWref_dq = dq_blk_size;
0310 matidx_dq = zeros(dq_blk_size,3);
0311 dq_ptr = 0;
0312 
0313 for i=1:nmaxclique
0314     nmaxcliquei = length(maxclique{i});
0315     for k=1:nmaxcliquei % row
0316         for m=k:nmaxcliquei % column
0317             
0318             % Check if this pair loop{i}(k) and loop{i}(m) has appeared in
0319             % any previous maxclique. If not, it isn't in Wref_dd.
0320             Wref_dd_found = 0;
0321             if i > 1
0322                 for r = 1:i-1
0323                    if any(maxclique{r}(:) == maxclique{i}(k)) && any(maxclique{r}(:) == maxclique{i}(m))
0324                        Wref_dd_found = 1;
0325                        break;
0326                    end
0327                 end
0328             end
0329 
0330             if ~Wref_dd_found
0331                 % If we didn't find this element already, add it to Wref_dd
0332                 % and matidx_dd
0333                 dd_ptr = dd_ptr + 1;
0334                 
0335                 if dd_ptr > lWref_dd
0336                     % Expand the variables by dd_blk_size
0337                     Wref_dd = [Wref_dd; zeros(dd_blk_size,2)];
0338                     matidx_dd = [matidx_dd; zeros(dd_ptr-1 + dd_blk_size,3)];
0339                     lWref_dd = length(Wref_dd(:,1));
0340                 end
0341                 
0342                 Wref_dd(dd_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0343                 matidx_dd(dd_ptr,1:3) = [i k m];
0344             end
0345 
0346             Wref_dq_found = Wref_dd_found;
0347 
0348             if ~Wref_dq_found
0349                 % If we didn't find this element already, add it to Wref_dq
0350                 % and matidx_dq
0351                 dq_ptr = dq_ptr + 1;
0352                 
0353                 if dq_ptr > lWref_dq
0354                     % Expand the variables by dq_blk_size
0355                     Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0356                     matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0357                     lWref_dq = length(Wref_dq(:,1));
0358                 end
0359                 
0360                 Wref_dq(dq_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0361                 matidx_dq(dq_ptr,1:3) = [i k m+nmaxcliquei];
0362             
0363                 if k ~= m % Already have the diagonal terms of the off-diagonal block
0364 
0365                     dq_ptr = dq_ptr + 1;
0366 
0367                     if dq_ptr > lWref_dq
0368                         % Expand the variables by dq_blk_size
0369                         Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0370                         matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0371                     end
0372 
0373                     Wref_dq(dq_ptr,1:2) = [maxclique{i}(m) maxclique{i}(k)];
0374                     matidx_dq(dq_ptr,1:3) = [i k+nmaxcliquei m];
0375                 end
0376             end
0377         end
0378     end
0379     
0380     if verbose >= 2 && mod(i,ndisplay) == 0
0381         fprintf('Loop identification: Loop %i of %i\n',i,nmaxclique);
0382     end
0383 end
0384 
0385 % Trim off excess zeros and empty matrix structures
0386 Wref_dd = Wref_dd(1:dd_ptr,:);
0387 matidx_dd = matidx_dd(1:dd_ptr,:);
0388 
0389 % Store index of Wref variables
0390 Wref_dd = [(1:length(Wref_dd)).' Wref_dd];
0391 Wref_dq = [(1:length(Wref_dq)).' Wref_dq];
0392 
0393 Wref_qq = Wref_dd;
0394 matidx_qq = zeros(size(Wref_qq,1),3);
0395 for i=1:size(Wref_qq,1) 
0396     nmaxcliquei = length(maxclique{matidx_dd(i,1)});
0397     matidx_qq(i,1:3) = [matidx_dd(i,1) matidx_dd(i,2)+nmaxcliquei matidx_dd(i,3)+nmaxcliquei];
0398 end
0399 
0400 
0401 %% Enforce linking constraints in dual
0402 
0403 for i=1:nmaxclique
0404     A{i} = [];
0405 end
0406 
0407 % Count the number of required linking constraints
0408 nbeta = 0;
0409 for i=1:size(E,1)
0410     overlap_idx = intersect(maxclique{E(i,1)},maxclique{E(i,2)});
0411     for k=1:length(overlap_idx)
0412         for m=k:length(overlap_idx)
0413             E11idx = find(maxclique{E(i,1)} == overlap_idx(k));
0414             E12idx = find(maxclique{E(i,1)} == overlap_idx(m));
0415             
0416             nbeta = nbeta + 3;
0417             
0418             if E11idx ~= E12idx
0419                 nbeta = nbeta + 1;
0420             end
0421         end
0422     end
0423     
0424     if verbose >= 2 && mod(i,ndisplay) == 0
0425         fprintf('Counting beta: %i of %i\n',i,size(E,1));
0426     end
0427 end
0428 
0429 % Make beta sdpvar
0430 if nbeta > 0
0431     beta = sdpvar(nbeta,1);
0432 end
0433 
0434 % Create the linking constraints defined by the maximal clique tree
0435 beta_idx = 0;
0436 for i=1:size(E,1)
0437     overlap_idx = intersect(maxclique{E(i,1)},maxclique{E(i,2)});
0438     nmaxclique1 = length(maxclique{E(i,1)});
0439     nmaxclique2 = length(maxclique{E(i,2)});
0440     for k=1:length(overlap_idx)
0441         for m=k:length(overlap_idx)
0442             E11idx = find(maxclique{E(i,1)} == overlap_idx(k));
0443             E12idx = find(maxclique{E(i,1)} == overlap_idx(m));
0444             E21idx = find(maxclique{E(i,2)} == overlap_idx(k));
0445             E22idx = find(maxclique{E(i,2)} == overlap_idx(m));
0446             
0447             beta_idx = beta_idx + 1;
0448             if ~isempty(A{E(i,1)})
0449                 A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0450             else
0451                 A{E(i,1)} = 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0452             end
0453             if ~isempty(A{E(i,2)})
0454                 A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0455             else
0456                 A{E(i,2)} = -0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0457             end
0458             
0459             beta_idx = beta_idx + 1;
0460             A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx+nmaxclique1; E12idx+nmaxclique1],[E12idx+nmaxclique1; E11idx+nmaxclique1], [1;1], 2*nmaxclique1, 2*nmaxclique1);
0461             A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx+nmaxclique2; E22idx+nmaxclique2],[E22idx+nmaxclique2; E21idx+nmaxclique2], [1;1], 2*nmaxclique2, 2*nmaxclique2);
0462             
0463             beta_idx = beta_idx + 1;
0464             A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx; E12idx+nmaxclique1], [E12idx+nmaxclique1; E11idx], [1;1], 2*nmaxclique1, 2*nmaxclique1);
0465             A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx; E22idx+nmaxclique2], [E22idx+nmaxclique2; E21idx], [1;1], 2*nmaxclique2, 2*nmaxclique2);
0466             
0467             if E11idx ~= E12idx
0468                 beta_idx = beta_idx + 1;
0469                 A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx+nmaxclique1; E12idx], [E12idx; E11idx+nmaxclique1], [1;1], 2*nmaxclique1, 2*nmaxclique1);
0470                 A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx+nmaxclique2; E22idx], [E22idx; E21idx+nmaxclique2], [1;1], 2*nmaxclique2, 2*nmaxclique2);
0471             end
0472         end
0473     end
0474     
0475     if verbose >= 2 && mod(i,ndisplay) == 0
0476         fprintf('SDP linking: %i of %i\n',i,size(E,1));
0477     end
0478     
0479 end
0480 
0481 % For systems with only one maximal clique, A doesn't get defined above.
0482 % Explicitly define it here.
0483 if nmaxclique == 1
0484     A{1} = sparse(2*nbus,2*nbus);
0485 end
0486 
0487 %% Build matrices
0488 
0489 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc);
0490 
0491 %% Create SDP variables for dual
0492 
0493 % lambda: active power equality constraints
0494 % gamma:  reactive power equality constraints
0495 % mu:     voltage magnitude constraints
0496 
0497 % Bus variables
0498 lambda = sdpvar(nbus,1);
0499 gamma = sdpvar(nbus,1);
0500 mu = sdpvar(nbus,1);
0501 
0502 
0503 
0504 %% Build bus constraints
0505 
0506 cost = 0;
0507 for k=1:nbus
0508         
0509     if mpc.bus(k,BUS_TYPE) == PQ
0510         % PQ bus has equality constraints on P and Q
0511         
0512         % Active power constraint
0513         A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0514         
0515         % Reactive power constraint
0516         A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -gamma(k), maxclique);
0517         
0518         % Cost function
0519         cost = cost + lambda(k)*Pinj(k) + gamma(k)*Qinj(k);
0520         
0521     elseif mpc.bus(k,BUS_TYPE) == PV
0522         % Scale the voltage at PV buses in constant proportion to the slack
0523         % bus voltage. Don't set the PV bus voltages to any specific value.
0524 
0525         % alpha is the square of the ratio between the voltage magnitudes
0526         % of the slack and PV bus, such that Vpv = alpha*Vslack
0527         alpha = Vmag(k)^2 / Vslack^2; 
0528         
0529         % Active power constraint
0530         A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0531         
0532         % Voltage magnitude constraint
0533         A = addToA(Mk(k)-alpha*Mk(slackbus_idx), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -mu(k), maxclique);
0534         
0535         % Cost function
0536         cost = cost + lambda(k)*Pinj(k);
0537         
0538     elseif mpc.bus(k,BUS_TYPE) == REF
0539         % Reference bus is unconstrained.
0540     else
0541         error('insolvablepf: Invalid bus type');
0542     end
0543    
0544     if verbose >= 2 && mod(k,ndisplay) == 0
0545         fprintf('SDP creation: Bus %i of %i\n',k,nbus);
0546     end
0547     
0548 end % Loop through all buses
0549 
0550 %% Incorporate objective function (minimize Vslack^2)
0551 
0552 A = addToA(Mk(slackbus_idx), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 1, maxclique);
0553 
0554 %% Formulate dual psd constraints
0555 
0556 for i=1:nmaxclique
0557     % Can multiply A by any non-zero scalar. This may affect numerics
0558     % of the solver.
0559     constraints = [constraints; 1*A{i} >= 0]; 
0560 end
0561 
0562 
0563 %% Solve the SDP
0564 
0565 % Preserve warning settings
0566 S = warning;
0567 
0568 % Run sdp solver
0569 sdpinfo = solvesdp(constraints, -cost, sdpopts); % Negative cost to convert maximization to minimization problem
0570 if sdpinfo.problem == 2 || sdpinfo.problem == -2 || sdpinfo.problem == -3
0571     error(yalmiperror(sdpinfo.problem));
0572 end
0573 if ~have_feature('octave') || have_feature('octave', 'vnum') >= 4.001
0574     %% (avoid bug in Octave 4.0.x, where warning state is left corrupted)
0575     warning(S);
0576 end
0577 
0578 if verbose >= 2
0579     fprintf('Solver exit message: %s\n',sdpinfo.info);
0580 end
0581 
0582 %% Calculate rank characteristics of the solution
0583 
0584 mineigratio = inf;
0585 for i=1:length(A)
0586     evl = eig(double(A{i}));
0587 
0588     % Calculate mineigratio
0589     eigA_ratio = abs(evl(3) / evl(1));
0590 
0591     if eigA_ratio < mineigratio
0592         mineigratio = eigA_ratio;
0593     end
0594 end
0595 
0596 
0597 %% Output results
0598 
0599 Vslack_min = sqrt(abs(double(cost)));
0600 insolvable = Vslack_min > Vmag(slackbus_idx);
0601 sigma = Vmag(slackbus_idx) / Vslack_min;
0602 eta = sigma^2;
insolvablepf

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
