0001 function [insolvable,Vslack_min,sigma,eta,mineigratio] = insolvablepf(mpc,mpopt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084 if nargin < 2
0085 mpopt = mpoption;
0086 end
0087
0088 mpc = loadcase(mpc);
0089 mpc = ext2int(mpc);
0090
0091
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;
0096 ndisplay = mpopt.sdp_pf.ndisplay;
0097 cholopts.dense = mpopt.sdp_pf.choldense;
0098 cholopts.aggressive = mpopt.sdp_pf.cholaggressive;
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
0109 sdpopts = yalmip_options([], mpopt);
0110
0111
0112
0113
0114
0115
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
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
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
0162
0163
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
0175
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
0191
0192
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
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
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
0245
0246
0247 if maxNumberOfCliques ~= 1 && nbus > 3
0248
0249
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
0255 [MC,ischordal] = maxCardSearch(Aadj);
0256
0257 if ~ischordal
0258
0259
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
0275
0276 if maxNumberOfCliques ~= 1 && nbus > 3
0277
0278
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
0288
0289
0290 end
0291 end
0292 cliqueCost = max(cliqueCost,cliqueCost.');
0293
0294
0295 cliqueCost = -cliqueCost;
0296 [E] = prim(cliqueCost);
0297
0298
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
0308
0309 tic;
0310
0311 constraints = [];
0312
0313
0314 dd_blk_size = 50*nbus;
0315 dq_blk_size = 2*dd_blk_size;
0316
0317 Wref_dd = zeros(dd_blk_size,2);
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);
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
0330 for m=k:nmaxcliquei
0331
0332
0333
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
0346
0347 dd_ptr = dd_ptr + 1;
0348
0349 if dd_ptr > lWref_dd
0350
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
0364
0365 dq_ptr = dq_ptr + 1;
0366
0367 if dq_ptr > lWref_dq
0368
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
0378
0379 dq_ptr = dq_ptr + 1;
0380
0381 if dq_ptr > lWref_dq
0382
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
0400 Wref_dd = Wref_dd(1:dd_ptr,:);
0401 matidx_dd = matidx_dd(1:dd_ptr,:);
0402
0403
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
0416
0417 for i=1:nmaxclique
0418 A{i} = [];
0419 end
0420
0421
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
0444 if nbeta > 0
0445 beta = sdpvar(nbeta,1);
0446 end
0447
0448
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
0496
0497 if nmaxclique == 1
0498 A{1} = sparse(2*nbus,2*nbus);
0499 end
0500
0501
0502
0503 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc);
0504
0505
0506
0507
0508
0509
0510
0511
0512 lambda = sdpvar(nbus,1);
0513 gamma = sdpvar(nbus,1);
0514 mu = sdpvar(nbus,1);
0515
0516
0517
0518
0519
0520 cost = 0;
0521 for k=1:nbus
0522
0523 if mpc.bus(k,BUS_TYPE) == PQ
0524
0525
0526
0527 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0528
0529
0530 A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -gamma(k), maxclique);
0531
0532
0533 cost = cost + lambda(k)*Pinj(k) + gamma(k)*Qinj(k);
0534
0535 elseif mpc.bus(k,BUS_TYPE) == PV
0536
0537
0538
0539
0540
0541 alpha = Vmag(k)^2 / Vslack^2;
0542
0543
0544 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0545
0546
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
0550 cost = cost + lambda(k)*Pinj(k);
0551
0552 elseif mpc.bus(k,BUS_TYPE) == REF
0553
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
0563
0564
0565
0566 A = addToA(Mk(slackbus_idx), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 1, maxclique);
0567
0568
0569
0570 for i=1:nmaxclique
0571
0572
0573 constraints = [constraints; 1*A{i} >= 0];
0574 end
0575
0576
0577
0578
0579
0580 S = warning;
0581
0582
0583 sdpinfo = solvesdp(constraints, -cost, sdpopts);
0584 warning(S);
0585
0586 if verbose >= 2
0587 fprintf('Solver exit message: %s\n',sdpinfo.info);
0588 end
0589
0590
0591
0592 mineigratio = inf;
0593 for i=1:length(A)
0594 evl = eig(double(A{i}));
0595
0596
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
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;