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 if nargin < 2
0066 mpopt = mpoption;
0067 end
0068
0069 mpc = loadcase(mpc);
0070 mpc = ext2int(mpc);
0071
0072
0073 ignore_angle_lim = mpopt.opf.ignore_angle_lim;
0074 verbose = mpopt.verbose;
0075 enforce_Qlimits = mpopt.pf.enforce_q_lims;
0076 maxNumberOfCliques = mpopt.sdp_pf.max_number_of_cliques;
0077 ndisplay = mpopt.sdp_pf.ndisplay;
0078 cholopts.dense = mpopt.sdp_pf.choldense;
0079 cholopts.aggressive = mpopt.sdp_pf.cholaggressive;
0080
0081 if enforce_Qlimits > 0
0082 enforce_Qlimits = 1;
0083 end
0084
0085 if ~have_fcn('yalmip')
0086 error('insolvablepf: The software package YALMIP is required to run insolvablepf. See http://users.isy.liu.se/johanl/yalmip/');
0087 end
0088
0089
0090 sdpopts = yalmip_options([], mpopt);
0091
0092
0093
0094
0095
0096
0097
0098 if enforce_Qlimits
0099 if verbose > 0
0100 fprintf('Generator reactive power limits are enforced. Using function insolvablepf_limitQ.\n');
0101 end
0102 [insolvable,eta,mineigratio] = insolvablepf_limitQ(mpc,mpopt);
0103 Vslack_min = nan;
0104 sigma = nan;
0105 return;
0106 end
0107
0108
0109
0110 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0111 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0112 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0113 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0114 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0115 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0116 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0117 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0118 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0119
0120
0121
0122 if isfield(mpc, 'userfcn') && length(mpc.userfcn) > 0 && ...
0123 isfield(mpc.userfcn, 'formulation')
0124 c = cellfun(@func2str, {mpc.userfcn.formulation.fcn}, 'UniformOutput', 0);
0125 if strfind(strcat(c{:}), 'userfcn_dcline_formulation')
0126 error('insolvablepf: DC lines are not implemented in insolvablepf');
0127 end
0128 end
0129
0130 if toggle_dcline(mpc, 'status')
0131 error('insolvablepf: DC lines are not implemented in insolvablepf');
0132 end
0133
0134 nbus = size(mpc.bus,1);
0135 ngen = size(mpc.gen,1);
0136 nbranch = size(mpc.branch,1);
0137
0138 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0139 warning('insolvablepf: Angle difference constraints are not implemented in SDP_PF. Ignoring angle difference constraints.');
0140 end
0141
0142
0143
0144
0145 for i=1:ngen
0146 busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0147 if isempty(busidx) || ~(mpc.bus(busidx,BUS_TYPE) == PV || mpc.bus(busidx,BUS_TYPE) == REF)
0148 mpc.bus(busidx,BUS_TYPE) = PV;
0149 if verbose >= 1
0150 warning('insolvablepf: Bus %s has generator(s) but is listed as a PQ bus. Changing to a PV bus.',int2str(busidx));
0151 end
0152 end
0153 end
0154
0155
0156
0157 for i=1:nbus
0158 if mpc.bus(i,BUS_TYPE) == PV
0159 genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0160 if isempty(genidx)
0161 mpc.bus(i,BUS_TYPE) = PQ;
0162 if verbose >= 1
0163 warning('insolvablepf: PV bus %i has no associated generator! Changing these buses to PQ.',i);
0164 end
0165 end
0166 end
0167 end
0168
0169
0170
0171
0172
0173
0174
0175 if maxNumberOfCliques ~= 1
0176 [Ainc] = makeIncidence(mpc.bus,mpc.branch);
0177 sparsity_pattern = abs(Ainc.'*Ainc) + eye(nbus,nbus);
0178 per = amd(sparsity_pattern,cholopts);
0179 [Rchol,p] = chol(sparsity_pattern(per,per),'lower');
0180
0181 if p ~= 0
0182 error('insolvablepf: sparsity_pattern not positive definite!');
0183 end
0184 else
0185 per = 1:nbus;
0186 end
0187
0188
0189 mpc.bus = mpc.bus(per,:);
0190 mpc.bus(:,BUS_I) = 1:nbus;
0191
0192 for i=1:ngen
0193 mpc.gen(i,GEN_BUS) = find(per == mpc.gen(i,GEN_BUS));
0194 end
0195
0196 [mpc.gen genidx] = sortrows(mpc.gen,1);
0197 mpc.gencost = mpc.gencost(genidx,:);
0198
0199 for i=1:nbranch
0200 mpc.branch(i,F_BUS) = find(per == mpc.branch(i,F_BUS));
0201 mpc.branch(i,T_BUS) = find(per == mpc.branch(i,T_BUS));
0202 end
0203
0204
0205 Sbase = mpc.baseMVA;
0206
0207
0208 Qinj = -mpc.bus(:,QD) / Sbase;
0209 Vmag = mpc.bus(:,VM);
0210
0211 Pd = mpc.bus(:,PD) / Sbase;
0212 Pg = zeros(nbus,1);
0213 for i=1:nbus
0214 genidx = find(mpc.gen(:,GEN_BUS) == i);
0215 if ~isempty(genidx)
0216 Pg(i) = sum(mpc.gen(genidx,PG)) / Sbase;
0217 Vmag(i) = mpc.gen(genidx(1),VG);
0218 end
0219 end
0220 Pinj = Pg - Pd;
0221
0222 slackbus_idx = find(mpc.bus(:,BUS_TYPE) == 3);
0223 Vslack = Vmag(slackbus_idx);
0224
0225
0226
0227
0228 if maxNumberOfCliques ~= 1 && nbus > 3
0229
0230
0231 [f,t] = find(Rchol - diag(diag(Rchol)));
0232 Aadj = sparse(f,t,ones(length(f),1),nbus,nbus);
0233 Aadj = max(Aadj,Aadj');
0234
0235
0236 [MC,ischordal] = maxCardSearch(Aadj);
0237
0238 if ~ischordal
0239
0240
0241 error('insolvablepf: Chordal extension adjacency matrix is not chordal!');
0242 end
0243
0244 for i=1:size(MC,2)
0245 maxclique{i} = find(MC(:,i));
0246 end
0247
0248 else
0249 maxclique{1} = 1:nbus;
0250 nmaxclique = 1;
0251 E = [];
0252 end
0253
0254
0255
0256
0257 if maxNumberOfCliques ~= 1 && nbus > 3
0258
0259
0260 nmaxclique = length(maxclique);
0261 cliqueCost = sparse(nmaxclique,nmaxclique);
0262 for i=1:nmaxclique
0263 maxcliquei = maxclique{i};
0264 for k=i+1:nmaxclique
0265
0266 cliqueCost(i,k) = sum(ismembc(maxcliquei,maxclique{k}));
0267
0268
0269
0270
0271 end
0272 end
0273 cliqueCost = max(cliqueCost,cliqueCost.');
0274
0275
0276 cliqueCost = -cliqueCost;
0277 [E] = prim(cliqueCost);
0278
0279
0280 if verbose >= 2
0281 [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,ndisplay);
0282 else
0283 [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,inf);
0284 end
0285 nmaxclique = length(maxclique);
0286 end
0287
0288
0289
0290 tic;
0291
0292 constraints = [];
0293
0294
0295 dd_blk_size = 50*nbus;
0296 dq_blk_size = 2*dd_blk_size;
0297
0298 Wref_dd = zeros(dd_blk_size,2);
0299 lWref_dd = dd_blk_size;
0300 matidx_dd = zeros(dd_blk_size,3);
0301 dd_ptr = 0;
0302
0303 Wref_dq = zeros(dq_blk_size,2);
0304 lWref_dq = dq_blk_size;
0305 matidx_dq = zeros(dq_blk_size,3);
0306 dq_ptr = 0;
0307
0308 for i=1:nmaxclique
0309 nmaxcliquei = length(maxclique{i});
0310 for k=1:nmaxcliquei
0311 for m=k:nmaxcliquei
0312
0313
0314
0315 Wref_dd_found = 0;
0316 if i > 1
0317 for r = 1:i-1
0318 if any(maxclique{r}(:) == maxclique{i}(k)) && any(maxclique{r}(:) == maxclique{i}(m))
0319 Wref_dd_found = 1;
0320 break;
0321 end
0322 end
0323 end
0324
0325 if ~Wref_dd_found
0326
0327
0328 dd_ptr = dd_ptr + 1;
0329
0330 if dd_ptr > lWref_dd
0331
0332 Wref_dd = [Wref_dd; zeros(dd_blk_size,2)];
0333 matidx_dd = [matidx_dd; zeros(dd_ptr-1 + dd_blk_size,3)];
0334 lWref_dd = length(Wref_dd(:,1));
0335 end
0336
0337 Wref_dd(dd_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0338 matidx_dd(dd_ptr,1:3) = [i k m];
0339 end
0340
0341 Wref_dq_found = Wref_dd_found;
0342
0343 if ~Wref_dq_found
0344
0345
0346 dq_ptr = dq_ptr + 1;
0347
0348 if dq_ptr > lWref_dq
0349
0350 Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0351 matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0352 lWref_dq = length(Wref_dq(:,1));
0353 end
0354
0355 Wref_dq(dq_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0356 matidx_dq(dq_ptr,1:3) = [i k m+nmaxcliquei];
0357
0358 if k ~= m
0359
0360 dq_ptr = dq_ptr + 1;
0361
0362 if dq_ptr > lWref_dq
0363
0364 Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0365 matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0366 end
0367
0368 Wref_dq(dq_ptr,1:2) = [maxclique{i}(m) maxclique{i}(k)];
0369 matidx_dq(dq_ptr,1:3) = [i k+nmaxcliquei m];
0370 end
0371 end
0372 end
0373 end
0374
0375 if verbose >= 2 && mod(i,ndisplay) == 0
0376 fprintf('Loop identification: Loop %i of %i\n',i,nmaxclique);
0377 end
0378 end
0379
0380
0381 Wref_dd = Wref_dd(1:dd_ptr,:);
0382 matidx_dd = matidx_dd(1:dd_ptr,:);
0383
0384
0385 Wref_dd = [(1:length(Wref_dd)).' Wref_dd];
0386 Wref_dq = [(1:length(Wref_dq)).' Wref_dq];
0387
0388 Wref_qq = Wref_dd;
0389 matidx_qq = zeros(size(Wref_qq,1),3);
0390 for i=1:size(Wref_qq,1)
0391 nmaxcliquei = length(maxclique{matidx_dd(i,1)});
0392 matidx_qq(i,1:3) = [matidx_dd(i,1) matidx_dd(i,2)+nmaxcliquei matidx_dd(i,3)+nmaxcliquei];
0393 end
0394
0395
0396
0397
0398 for i=1:nmaxclique
0399 A{i} = [];
0400 end
0401
0402
0403 nbeta = 0;
0404 for i=1:size(E,1)
0405 overlap_idx = intersect(maxclique{E(i,1)},maxclique{E(i,2)});
0406 for k=1:length(overlap_idx)
0407 for m=k:length(overlap_idx)
0408 E11idx = find(maxclique{E(i,1)} == overlap_idx(k));
0409 E12idx = find(maxclique{E(i,1)} == overlap_idx(m));
0410
0411 nbeta = nbeta + 3;
0412
0413 if E11idx ~= E12idx
0414 nbeta = nbeta + 1;
0415 end
0416 end
0417 end
0418
0419 if verbose >= 2 && mod(i,ndisplay) == 0
0420 fprintf('Counting beta: %i of %i\n',i,size(E,1));
0421 end
0422 end
0423
0424
0425 if nbeta > 0
0426 beta = sdpvar(nbeta,1);
0427 end
0428
0429
0430 beta_idx = 0;
0431 for i=1:size(E,1)
0432 overlap_idx = intersect(maxclique{E(i,1)},maxclique{E(i,2)});
0433 nmaxclique1 = length(maxclique{E(i,1)});
0434 nmaxclique2 = length(maxclique{E(i,2)});
0435 for k=1:length(overlap_idx)
0436 for m=k:length(overlap_idx)
0437 E11idx = find(maxclique{E(i,1)} == overlap_idx(k));
0438 E12idx = find(maxclique{E(i,1)} == overlap_idx(m));
0439 E21idx = find(maxclique{E(i,2)} == overlap_idx(k));
0440 E22idx = find(maxclique{E(i,2)} == overlap_idx(m));
0441
0442 beta_idx = beta_idx + 1;
0443 if ~isempty(A{E(i,1)})
0444 A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0445 else
0446 A{E(i,1)} = 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0447 end
0448 if ~isempty(A{E(i,2)})
0449 A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0450 else
0451 A{E(i,2)} = -0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0452 end
0453
0454 beta_idx = beta_idx + 1;
0455 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);
0456 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);
0457
0458 beta_idx = beta_idx + 1;
0459 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);
0460 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);
0461
0462 if E11idx ~= E12idx
0463 beta_idx = beta_idx + 1;
0464 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);
0465 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);
0466 end
0467 end
0468 end
0469
0470 if verbose >= 2 && mod(i,ndisplay) == 0
0471 fprintf('SDP linking: %i of %i\n',i,size(E,1));
0472 end
0473
0474 end
0475
0476
0477
0478 if nmaxclique == 1
0479 A{1} = sparse(2*nbus,2*nbus);
0480 end
0481
0482
0483
0484 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc);
0485
0486
0487
0488
0489
0490
0491
0492
0493 lambda = sdpvar(nbus,1);
0494 gamma = sdpvar(nbus,1);
0495 mu = sdpvar(nbus,1);
0496
0497
0498
0499
0500
0501 cost = 0;
0502 for k=1:nbus
0503
0504 if mpc.bus(k,BUS_TYPE) == PQ
0505
0506
0507
0508 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0509
0510
0511 A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -gamma(k), maxclique);
0512
0513
0514 cost = cost + lambda(k)*Pinj(k) + gamma(k)*Qinj(k);
0515
0516 elseif mpc.bus(k,BUS_TYPE) == PV
0517
0518
0519
0520
0521
0522 alpha = Vmag(k)^2 / Vslack^2;
0523
0524
0525 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambda(k), maxclique);
0526
0527
0528 A = addToA(Mk(k)-alpha*Mk(slackbus_idx), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -mu(k), maxclique);
0529
0530
0531 cost = cost + lambda(k)*Pinj(k);
0532
0533 elseif mpc.bus(k,BUS_TYPE) == REF
0534
0535 else
0536 error('insolvablepf: Invalid bus type');
0537 end
0538
0539 if verbose >= 2 && mod(k,ndisplay) == 0
0540 fprintf('SDP creation: Bus %i of %i\n',k,nbus);
0541 end
0542
0543 end
0544
0545
0546
0547 A = addToA(Mk(slackbus_idx), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 1, maxclique);
0548
0549
0550
0551 for i=1:nmaxclique
0552
0553
0554 constraints = [constraints; 1*A{i} >= 0];
0555 end
0556
0557
0558
0559
0560
0561 S = warning;
0562
0563
0564 sdpinfo = solvesdp(constraints, -cost, sdpopts);
0565 warning(S);
0566
0567 if verbose >= 2
0568 fprintf('Solver exit message: %s\n',sdpinfo.info);
0569 end
0570
0571
0572
0573 mineigratio = inf;
0574 for i=1:length(A)
0575 evl = eig(double(A{i}));
0576
0577
0578 eigA_ratio = abs(evl(3) / evl(1));
0579
0580 if eigA_ratio < mineigratio
0581 mineigratio = eigA_ratio;
0582 end
0583 end
0584
0585
0586
0587
0588 Vslack_min = sqrt(abs(double(cost)));
0589 insolvable = Vslack_min > Vmag(slackbus_idx);
0590 sigma = Vmag(slackbus_idx) / Vslack_min;
0591 eta = sigma^2;