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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
