Home > matpower5.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 %   $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;

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005