Home > matpower6.0 > extras > sdp_pf > insolvablepf.m

insolvablepf

PURPOSE ^

INSOLVABLEPF A sufficient condition for power flow insolvability

SYNOPSIS ^

function [insolvable,Vslack_min,sigma,eta,mineigratio] = insolvablepf(mpc,mpopt)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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;

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005