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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
