insolvablepf

INSOLVABLEPF A sufficient condition for power flow insolvability

   [INSOLVABLE,VSLACK_MIN,SIGMA,ETA,MINEIGRATIO] = INSOLVABLEPF(MPC,MPOPT)

   Evaluates a sufficient condition for insolvability of the power flow
   equations. The voltage at the slack bus is minimized (with 
   proportional changes in voltage magnitudes at PV buses) using a 
   semidefinite programming relaxation of the power flow equations. If the
   minimum achievable slack bus voltage is greater than the specified slack
   bus voltage, no power flow solutions exist. The converse does not
   necessarily hold; a power flow solution may not exist for cases
   where the output insolvable is equal to 0. See [1] and [2] for further
   details.

   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).
       VSLACK_MIN : Minimum possible slack voltage obtained from the
           semidefinite programming relaxation. The power flow equations
           are insolvable if Vslack_min > V0, where V0 is the specified
           voltage at the slack bus.
       SIGMA : Controlled voltage margin to the power flow solvability
           boundary.
       ETA : Power injection margin to the power flow solvability
           boundary in the profile of a uniform, constant power factor
           change in power injections.
       MINEIGRATIO : A measure of satisfaction of the rank relaxation.
           Large values indicate satisfaction. (Note that satisfaction of
           the rank relaxation is not required for correctness of the
           insolvability condition).

   Note that this function uses a matrix completion decomposition and is
   therefore suitable for large systems. 

 [1] D.K. Molzahn, B.C. Lesieutre, and C.L. DeMarco, "A Sufficient Condition
     for Power Flow Insolvability with Applications to Voltage Stability
     Margins," IEEE Transactions on Power Systems, vol. 28, no. 3,
     pp. 2592-2601, August 2013.

 [2] D.K. Molzahn, B.C. Lesieutre, and C.L. DeMarco, "A Sufficient Condition 
     for Power Flow Insolvability with Applications to Voltage Stability 
     Margins," University of Wisconsin-Madison Department of Electrical
     and Computer Engineering, Tech. Rep. ECE-12-01, 2012, [Online]. 
     Available: http://arxiv.org/abs/1204.6285.

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