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