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 if nargin < 2
0064 mpopt = mpoption;
0065 end
0066
0067
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
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;
0082 ndisplay = mpopt.sdp_pf.ndisplay;
0083 cholopts.dense = mpopt.sdp_pf.choldense;
0084 cholopts.aggressive = mpopt.sdp_pf.cholaggressive;
0085
0086 if enforce_Qlimits > 0
0087 enforce_Qlimits = 1;
0088 end
0089
0090 if ~have_feature('yalmip')
0091 error('insolvablepf: The software package YALMIP is required to run insolvablepf. See https://yalmip.github.io');
0092 end
0093
0094
0095 sdpopts = yalmip_options([], mpopt);
0096
0097
0098
0099
0100
0101
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
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
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
0148
0149
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
0161
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
0177
0178
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
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
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
0231
0232
0233 if maxNumberOfCliques ~= 1 && nbus > 3
0234
0235
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
0241 [MC,ischordal] = maxCardSearch(Aadj);
0242
0243 if ~ischordal
0244
0245
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
0261
0262 if maxNumberOfCliques ~= 1 && nbus > 3
0263
0264
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
0274
0275
0276 end
0277 end
0278 cliqueCost = max(cliqueCost,cliqueCost.');
0279
0280
0281 cliqueCost = -cliqueCost;
0282 [E] = prim(cliqueCost);
0283
0284
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
0294
0295 tic;
0296
0297 constraints = [];
0298
0299
0300 dd_blk_size = 50*nbus;
0301 dq_blk_size = 2*dd_blk_size;
0302
0303 Wref_dd = zeros(dd_blk_size,2);
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);
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
0316 for m=k:nmaxcliquei
0317
0318
0319
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
0332
0333 dd_ptr = dd_ptr + 1;
0334
0335 if dd_ptr > lWref_dd
0336
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
0350
0351 dq_ptr = dq_ptr + 1;
0352
0353 if dq_ptr > lWref_dq
0354
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
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 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
0386 Wref_dd = Wref_dd(1:dd_ptr,:);
0387 matidx_dd = matidx_dd(1:dd_ptr,:);
0388
0389
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
0402
0403 for i=1:nmaxclique
0404 A{i} = [];
0405 end
0406
0407
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
0430 if nbeta > 0
0431 beta = sdpvar(nbeta,1);
0432 end
0433
0434
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
0482
0483 if nmaxclique == 1
0484 A{1} = sparse(2*nbus,2*nbus);
0485 end
0486
0487
0488
0489 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc);
0490
0491
0492
0493
0494
0495
0496
0497
0498 lambda = sdpvar(nbus,1);
0499 gamma = sdpvar(nbus,1);
0500 mu = sdpvar(nbus,1);
0501
0502
0503
0504
0505
0506 cost = 0;
0507 for k=1:nbus
0508
0509 if mpc.bus(k,BUS_TYPE) == PQ
0510
0511
0512
0513 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0514
0515
0516 A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -gamma(k), maxclique);
0517
0518
0519 cost = cost + lambda(k)*Pinj(k) + gamma(k)*Qinj(k);
0520
0521 elseif mpc.bus(k,BUS_TYPE) == PV
0522
0523
0524
0525
0526
0527 alpha = Vmag(k)^2 / Vslack^2;
0528
0529
0530 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0531
0532
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
0536 cost = cost + lambda(k)*Pinj(k);
0537
0538 elseif mpc.bus(k,BUS_TYPE) == REF
0539
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
0549
0550
0551
0552 A = addToA(Mk(slackbus_idx), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 1, maxclique);
0553
0554
0555
0556 for i=1:nmaxclique
0557
0558
0559 constraints = [constraints; 1*A{i} >= 0];
0560 end
0561
0562
0563
0564
0565
0566 S = warning;
0567
0568
0569 sdpinfo = solvesdp(constraints, -cost, sdpopts);
0570 if sdpinfo.problem == 2 || sdpinfo.problem == -2 || sdpinfo.problem == -3
0571 error(yalmiperror(sdpinfo.problem));
0572 end
0573 if ~have_feature('octave') || have_feature('octave', 'vnum') >= 4.001
0574
0575 warning(S);
0576 end
0577
0578 if verbose >= 2
0579 fprintf('Solver exit message: %s\n',sdpinfo.info);
0580 end
0581
0582
0583
0584 mineigratio = inf;
0585 for i=1:length(A)
0586 evl = eig(double(A{i}));
0587
0588
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
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;