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

Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005