0001 function [results,success,raw] = sdpopf_solver(om, 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 ignore_angle_lim = mpopt.opf.ignore_angle_lim;
0049 verbose = mpopt.verbose;
0050 maxNumberOfCliques = mpopt.sdp_pf.max_number_of_cliques;
0051 eps_r = mpopt.sdp_pf.eps_r;
0052 recover_voltage = mpopt.sdp_pf.recover_voltage;
0053 recover_injections = mpopt.sdp_pf.recover_injections;
0054 minPgendiff = mpopt.sdp_pf.min_Pgen_diff;
0055 minQgendiff = mpopt.sdp_pf.min_Qgen_diff;
0056 maxlinelimit = mpopt.sdp_pf.max_line_limit;
0057 maxgenlimit = mpopt.sdp_pf.max_gen_limit;
0058 ndisplay = mpopt.sdp_pf.ndisplay;
0059 cholopts.dense = mpopt.sdp_pf.choldense;
0060 cholopts.aggressive = mpopt.sdp_pf.cholaggressive;
0061 binding_lagrange = mpopt.sdp_pf.bind_lagrange;
0062 zeroeval_tol = mpopt.sdp_pf.zeroeval_tol;
0063 mineigratio_tol = mpopt.sdp_pf.mineigratio_tol;
0064
0065 if upper(mpopt.opf.flow_lim(1)) ~= 'S' && upper(mpopt.opf.flow_lim(1)) ~= 'P'
0066 error('sdpopf_solver: Only ''S'' and ''P'' options are currently implemented for MPOPT.opf.flow_lim when MPOPT.opf.ac.solver == ''SDPOPF''');
0067 end
0068
0069
0070 sdpopts = yalmip_options([], mpopt);
0071
0072 if verbose > 0
0073 v = sdp_pf_ver('all');
0074 fprintf('SDPOPF Version %s, %s', v.Version, v.Date);
0075 if ~isempty(sdpopts.solver)
0076 fprintf(' -- Using %s.\n\n',upper(sdpopts.solver));
0077 else
0078 fprintf('\n\n');
0079 end
0080 end
0081
0082 tic;
0083
0084
0085 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0086 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0087 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0088 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0089 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0090 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0091 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0092 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0093 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0094
0095
0096
0097
0098 mpc = get_mpc(om);
0099
0100 if toggle_dcline(mpc, 'status')
0101 error('sdpopf_solver: DC lines are not implemented in SDP_PF');
0102 end
0103
0104 nbus = size(mpc.bus,1);
0105 ngen = size(mpc.gen,1);
0106 nbranch = size(mpc.branch,1);
0107 bustype = mpc.bus(:,BUS_TYPE);
0108
0109 if any(mpc.gencost(:,NCOST) > 3 & mpc.gencost(:,MODEL) == POLYNOMIAL)
0110 error('sdpopf_solver: SDPOPF is limited to quadratic cost functions.');
0111 end
0112
0113 if size(mpc.gencost,1) > ngen && verbose >= 1
0114 warning('sdpopf_solver: Reactive power costs are not implemented in SDPOPF. Ignoring reactive power costs.');
0115 end
0116
0117 if ~ignore_angle_lim && (any(mpc.branch(:,ANGMIN) ~= -360) || any(mpc.branch(:,ANGMAX) ~= 360))
0118 warning('sdpopf_solver: Angle difference constraints are not implemented in SDPOPF. Ignoring angle difference constraints.');
0119 end
0120
0121
0122 mpc.branch(:,BR_R) = max(mpc.branch(:,BR_R),eps_r);
0123
0124
0125
0126
0127 for i=1:ngen
0128 busidx = find(mpc.bus(:,BUS_I) == mpc.gen(i,GEN_BUS));
0129 if isempty(busidx) || ~(bustype(busidx) == PV || bustype(busidx) == REF)
0130 bustype(busidx) = PV;
0131 end
0132 end
0133
0134
0135
0136 for i=1:nbus
0137 if bustype(i) == PV
0138 genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I), 1);
0139 if isempty(genidx)
0140 bustype(i) = PQ;
0141 end
0142 end
0143 end
0144
0145
0146
0147
0148
0149
0150
0151 for i=1:nbus
0152 genidx = find(mpc.gen(:,GEN_BUS) == mpc.bus(i,BUS_I));
0153 if ~all(mpc.gencost(genidx,MODEL) == PW_LINEAR) && ~all(mpc.gencost(genidx,MODEL) == POLYNOMIAL)
0154 error('sdpopf_solver: Bus %i has generators with both piecewise-linear and quadratic generators, which is not supported in this version of SDPOPF.',i);
0155 end
0156 end
0157
0158
0159 mpc.bus(:,size(mpc.bus,2)+1) = mpc.bus(:,BUS_I);
0160
0161 tsetup = toc;
0162
0163 tic;
0164
0165
0166
0167
0168 if maxNumberOfCliques ~= 1
0169 [Ainc] = makeIncidence(mpc.bus,mpc.branch);
0170 sparsity_pattern = abs(Ainc.'*Ainc) + speye(nbus,nbus);
0171 per = amd(sparsity_pattern,cholopts);
0172 [Rchol,p] = chol(sparsity_pattern(per,per),'lower');
0173
0174 if p ~= 0
0175 error('sdpopf_solver: sparsity_pattern not positive definite!');
0176 end
0177 else
0178 per = 1:nbus;
0179 end
0180
0181
0182 mpc.bus = mpc.bus(per,:);
0183 mpc.bus(:,BUS_I) = 1:nbus;
0184 bustype = bustype(per);
0185
0186 for i=1:ngen
0187 mpc.gen(i,GEN_BUS) = find(per == mpc.gen(i,GEN_BUS));
0188 end
0189
0190 [mpc.gen genidx] = sortrows(mpc.gen,1);
0191 mpc.gencost = mpc.gencost(genidx,:);
0192
0193 for i=1:nbranch
0194 mpc.branch(i,F_BUS) = find(per == mpc.branch(i,F_BUS));
0195 mpc.branch(i,T_BUS) = find(per == mpc.branch(i,T_BUS));
0196 end
0197
0198
0199 Sbase = mpc.baseMVA;
0200
0201 Pd = mpc.bus(:,PD) / Sbase;
0202 Qd = mpc.bus(:,QD) / Sbase;
0203
0204 for i=1:nbus
0205 genidx = find(mpc.gen(:,GEN_BUS) == i);
0206 if ~isempty(genidx)
0207 Qmax(i) = sum(mpc.gen(genidx,QMAX)) / Sbase;
0208 Qmin(i) = sum(mpc.gen(genidx,QMIN)) / Sbase;
0209 else
0210 Qmax(i) = 0;
0211 Qmin(i) = 0;
0212 end
0213 end
0214
0215 Vmax = mpc.bus(:,VMAX);
0216 Vmin = mpc.bus(:,VMIN);
0217
0218 Smax = mpc.branch(:,RATE_A) / Sbase;
0219
0220
0221
0222 c2 = zeros(ngen,1);
0223 c1 = zeros(ngen,1);
0224 c0 = zeros(ngen,1);
0225 for genidx=1:ngen
0226 if mpc.gencost(genidx,MODEL) == POLYNOMIAL
0227 if mpc.gencost(genidx,NCOST) == 3
0228 c2(genidx) = mpc.gencost(genidx,COST) * Sbase^2;
0229 c1(genidx) = mpc.gencost(genidx,COST+1) * Sbase^1;
0230 c0(genidx) = mpc.gencost(genidx,COST+2) * Sbase^0;
0231 elseif mpc.gencost(genidx,NCOST) == 2
0232 c1(genidx) = mpc.gencost(genidx,COST) * Sbase^1;
0233 c0(genidx) = mpc.gencost(genidx,COST+1) * Sbase^0;
0234 else
0235 error('sdpopf_solver: Only piecewise-linear and quadratic cost functions are implemented in SDPOPF');
0236 end
0237 end
0238 end
0239
0240 if any(c2 < 0)
0241 error('sdpopf_solver: Quadratic term of generator cost function is negative. Must be convex (non-negative).');
0242 end
0243
0244 maxlinelimit = maxlinelimit / Sbase;
0245 maxgenlimit = maxgenlimit / Sbase;
0246 minPgendiff = minPgendiff / Sbase;
0247 minQgendiff = minQgendiff / Sbase;
0248
0249
0250 [Y, Yf, Yt] = makeYbus(mpc);
0251
0252
0253
0254
0255
0256 if maxNumberOfCliques ~= 1 && nbus > 3
0257
0258
0259 [f,t] = find(Rchol - diag(diag(Rchol)));
0260 Aadj = sparse(f,t,ones(length(f),1),nbus,nbus);
0261 Aadj = max(Aadj,Aadj');
0262
0263
0264 [MC,ischordal] = maxCardSearch(Aadj);
0265
0266 if ~ischordal
0267
0268
0269 error('sdpopf_solver: Chordal extension adjacency matrix is not chordal!');
0270 end
0271
0272 for i=1:size(MC,2)
0273 maxclique{i} = find(MC(:,i));
0274 end
0275
0276 else
0277 maxclique{1} = 1:nbus;
0278 nmaxclique = 1;
0279 E = [];
0280 end
0281
0282
0283
0284
0285 if maxNumberOfCliques ~= 1 && nbus > 3
0286
0287
0288 nmaxclique = length(maxclique);
0289 cliqueCost = sparse(nmaxclique,nmaxclique);
0290 for i=1:nmaxclique
0291 maxcliquei = maxclique{i};
0292 for k=i+1:nmaxclique
0293
0294 cliqueCost(i,k) = sum(ismembc(maxcliquei,maxclique{k}));
0295
0296
0297
0298
0299 end
0300 end
0301 cliqueCost = max(cliqueCost,cliqueCost.');
0302
0303
0304 cliqueCost = -cliqueCost;
0305 [E] = prim(cliqueCost);
0306
0307
0308 if verbose >= 2
0309 [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,ndisplay);
0310 else
0311 [maxclique,E] = combineMaxCliques(maxclique,E,maxNumberOfCliques,inf);
0312 end
0313 nmaxclique = length(maxclique);
0314 end
0315 tmaxclique = toc;
0316
0317
0318
0319 tic;
0320
0321 constraints = [];
0322
0323
0324 dd_blk_size = 50*nbus;
0325 dq_blk_size = 2*dd_blk_size;
0326
0327 Wref_dd = zeros(dd_blk_size,2);
0328 lWref_dd = dd_blk_size;
0329 matidx_dd = zeros(dd_blk_size,3);
0330 dd_ptr = 0;
0331
0332 Wref_dq = zeros(dq_blk_size,2);
0333 lWref_dq = dq_blk_size;
0334 matidx_dq = zeros(dq_blk_size,3);
0335 dq_ptr = 0;
0336
0337 for i=1:nmaxclique
0338 nmaxcliquei = length(maxclique{i});
0339 for k=1:nmaxcliquei
0340 for m=k:nmaxcliquei
0341
0342
0343
0344 Wref_dd_found = 0;
0345 if i > 1
0346 for r = 1:i-1
0347 if any(maxclique{r}(:) == maxclique{i}(k)) && any(maxclique{r}(:) == maxclique{i}(m))
0348 Wref_dd_found = 1;
0349 break;
0350 end
0351 end
0352 end
0353
0354 if ~Wref_dd_found
0355
0356
0357 dd_ptr = dd_ptr + 1;
0358
0359 if dd_ptr > lWref_dd
0360
0361 Wref_dd = [Wref_dd; zeros(dd_blk_size,2)];
0362 matidx_dd = [matidx_dd; zeros(dd_ptr-1 + dd_blk_size,3)];
0363 lWref_dd = length(Wref_dd(:,1));
0364 end
0365
0366 Wref_dd(dd_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0367 matidx_dd(dd_ptr,1:3) = [i k m];
0368 end
0369
0370 Wref_dq_found = Wref_dd_found;
0371
0372 if ~Wref_dq_found
0373
0374
0375 dq_ptr = dq_ptr + 1;
0376
0377 if dq_ptr > lWref_dq
0378
0379 Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0380 matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0381 lWref_dq = length(Wref_dq(:,1));
0382 end
0383
0384 Wref_dq(dq_ptr,1:2) = [maxclique{i}(k) maxclique{i}(m)];
0385 matidx_dq(dq_ptr,1:3) = [i k m+nmaxcliquei];
0386
0387 if k ~= m
0388
0389 dq_ptr = dq_ptr + 1;
0390
0391 if dq_ptr > lWref_dq
0392
0393 Wref_dq = [Wref_dq; zeros(dq_blk_size,2)];
0394 matidx_dq = [matidx_dq; zeros(dq_ptr-1 + dq_blk_size,3)];
0395 end
0396
0397 Wref_dq(dq_ptr,1:2) = [maxclique{i}(m) maxclique{i}(k)];
0398 matidx_dq(dq_ptr,1:3) = [i k+nmaxcliquei m];
0399 end
0400 end
0401 end
0402 end
0403
0404 if verbose >= 2 && mod(i,ndisplay) == 0
0405 fprintf('Loop identification: Loop %i of %i\n',i,nmaxclique);
0406 end
0407 end
0408
0409
0410 Wref_dd = Wref_dd(1:dd_ptr,:);
0411 matidx_dd = matidx_dd(1:dd_ptr,:);
0412
0413
0414 Wref_dd = [(1:length(Wref_dd)).' Wref_dd];
0415 Wref_dq = [(1:length(Wref_dq)).' Wref_dq];
0416
0417 Wref_qq = Wref_dd;
0418 matidx_qq = zeros(size(Wref_qq,1),3);
0419 for i=1:size(Wref_qq,1)
0420 nmaxcliquei = length(maxclique{matidx_dd(i,1)});
0421 matidx_qq(i,1:3) = [matidx_dd(i,1) matidx_dd(i,2)+nmaxcliquei matidx_dd(i,3)+nmaxcliquei];
0422 end
0423
0424
0425
0426
0427 for i=1:nmaxclique
0428 A{i} = [];
0429 end
0430
0431
0432 nbeta = 0;
0433 for i=1:size(E,1)
0434 overlap_idx = intersect(maxclique{E(i,1)},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
0440 nbeta = nbeta + 3;
0441
0442 if E11idx ~= E12idx
0443 nbeta = nbeta + 1;
0444 end
0445 end
0446 end
0447
0448 if verbose >= 2 && mod(i,ndisplay) == 0
0449 fprintf('Counting beta: %i of %i\n',i,size(E,1));
0450 end
0451 end
0452
0453
0454 if nbeta > 0
0455 beta = sdpvar(nbeta,1);
0456 end
0457
0458
0459 beta_idx = 0;
0460 for i=1:size(E,1)
0461 overlap_idx = intersect(maxclique{E(i,1)},maxclique{E(i,2)});
0462 nmaxclique1 = length(maxclique{E(i,1)});
0463 nmaxclique2 = length(maxclique{E(i,2)});
0464 for k=1:length(overlap_idx)
0465 for m=k:length(overlap_idx)
0466 E11idx = find(maxclique{E(i,1)} == overlap_idx(k));
0467 E12idx = find(maxclique{E(i,1)} == overlap_idx(m));
0468 E21idx = find(maxclique{E(i,2)} == overlap_idx(k));
0469 E22idx = find(maxclique{E(i,2)} == overlap_idx(m));
0470
0471 beta_idx = beta_idx + 1;
0472 if ~isempty(A{E(i,1)})
0473 A{E(i,1)} = A{E(i,1)} + 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0474 else
0475 A{E(i,1)} = 0.5*beta(beta_idx)*sparse([E11idx; E12idx], [E12idx; E11idx], [1; 1], 2*nmaxclique1, 2*nmaxclique1);
0476 end
0477 if ~isempty(A{E(i,2)})
0478 A{E(i,2)} = A{E(i,2)} - 0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0479 else
0480 A{E(i,2)} = -0.5*beta(beta_idx)*sparse([E21idx; E22idx], [E22idx; E21idx], [1; 1], 2*nmaxclique2, 2*nmaxclique2);
0481 end
0482
0483 beta_idx = beta_idx + 1;
0484 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);
0485 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);
0486
0487 beta_idx = beta_idx + 1;
0488 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);
0489 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);
0490
0491 if E11idx ~= E12idx
0492 beta_idx = beta_idx + 1;
0493 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);
0494 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);
0495 end
0496 end
0497 end
0498
0499 if verbose >= 2 && mod(i,ndisplay) == 0
0500 fprintf('SDP linking: %i of %i\n',i,size(E,1));
0501 end
0502
0503 end
0504
0505
0506
0507 if nmaxclique == 1
0508 A{1} = sparse(2*nbus,2*nbus);
0509 end
0510
0511
0512
0513 [Yk,Yk_,Mk,Ylineft,Ylinetf,Y_lineft,Y_linetf] = makesdpmat(mpc);
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523 nlambda_eq = 0;
0524 ngamma_eq = 0;
0525 ngamma_ineq = 0;
0526 nmu_ineq = 0;
0527
0528 for i=1:nbus
0529
0530 nlambda_eq = nlambda_eq + 1;
0531 lambda_eq(i) = nlambda_eq;
0532
0533 if bustype(i) == PQ || (Qmax(i) - Qmin(i) < minQgendiff)
0534 ngamma_eq = ngamma_eq + 1;
0535 gamma_eq(i) = ngamma_eq;
0536 else
0537 gamma_eq(i) = nan;
0538 end
0539
0540 if (bustype(i) == PV || bustype(i) == REF) && (Qmax(i) - Qmin(i) >= minQgendiff) && (Qmax(i) <= maxgenlimit || Qmin(i) >= -maxgenlimit)
0541 ngamma_ineq = ngamma_ineq + 1;
0542 gamma_ineq(i) = ngamma_ineq;
0543 else
0544 gamma_ineq(i) = nan;
0545 end
0546
0547 nmu_ineq = nmu_ineq + 1;
0548 mu_ineq(i) = nmu_ineq;
0549
0550 psiU_sdp{i} = [];
0551 psiL_sdp{i} = [];
0552
0553 pwl_sdp{i} = [];
0554 end
0555
0556 lambdaeq_sdp = sdpvar(nlambda_eq,1);
0557
0558 gammaeq_sdp = sdpvar(ngamma_eq,1);
0559 gammaU_sdp = sdpvar(ngamma_ineq,1);
0560 gammaL_sdp = sdpvar(ngamma_ineq,1);
0561
0562 muU_sdp = sdpvar(nmu_ineq,1);
0563 muL_sdp = sdpvar(nmu_ineq,1);
0564
0565
0566
0567 if upper(mpopt.opf.flow_lim(1)) == 'S'
0568 Hlookup = zeros(2*sum(Smax ~= 0 & Smax < maxlinelimit),2);
0569 end
0570
0571 nlconstraint = 0;
0572 for i=1:nbranch
0573 if Smax(i) ~= 0 && Smax(i) < maxlinelimit
0574 nlconstraint = nlconstraint + 1;
0575 Hlookup(nlconstraint,:) = [i 1];
0576 nlconstraint = nlconstraint + 1;
0577 Hlookup(nlconstraint,:) = [i 0];
0578 end
0579 end
0580
0581 if upper(mpopt.opf.flow_lim(1)) == 'S'
0582 for i=1:nlconstraint
0583 Hsdp{i} = sdpvar(3,3);
0584 end
0585 elseif upper(mpopt.opf.flow_lim(1)) == 'P'
0586 Hsdp = sdpvar(2*nlconstraint,1);
0587 end
0588
0589
0590
0591
0592
0593 for k=1:nbus
0594
0595 if bustype(k) == PQ
0596
0597
0598
0599
0600 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -lambdaeq_sdp(lambda_eq(k),1), maxclique);
0601
0602
0603
0604
0605
0606
0607
0608 A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, -gammaeq_sdp(gamma_eq(k),1), maxclique);
0609
0610
0611
0612
0613
0614
0615
0616 if exist('cost','var')
0617 cost = cost ...
0618 + lambdaeq_sdp(lambda_eq(k))*(-Pd(k)) ...
0619 + gammaeq_sdp(gamma_eq(k))*(-Qd(k));
0620 else
0621 cost = lambdaeq_sdp(lambda_eq(k))*(-Pd(k)) ...
0622 + gammaeq_sdp(gamma_eq(k))*(-Qd(k));
0623 end
0624
0625
0626 elseif bustype(k) == PV || bustype(k) == REF
0627
0628
0629
0630
0631
0632 genidx = find(mpc.gen(:,GEN_BUS) == k);
0633
0634 if isempty(genidx)
0635 error('sdpopf_solver: Generator for bus %i not found!',k);
0636 end
0637
0638
0639
0640
0641 psiU_sdp{k} = sdpvar(length(genidx),1);
0642 psiL_sdp{k} = sdpvar(length(genidx),1);
0643
0644 clear lambdai Pmax Pmin
0645 for i=1:length(genidx)
0646
0647
0648 if mpc.gencost(genidx(i),MODEL) == PW_LINEAR
0649 pwl = 1;
0650
0651 pwl_sdp{k}{i} = sdpvar(mpc.gencost(genidx(i),NCOST)-1,1);
0652
0653
0654 x_pw = mpc.gencost(genidx(i),-1+COST+(1:2:2*length(pwl_sdp{k}{i})+1)) / Sbase;
0655 c_pw = mpc.gencost(genidx(i),-1+COST+(2:2:2*length(pwl_sdp{k}{i})+2));
0656 else
0657 pwl = 0;
0658 pwl_sdp{k}{i} = [];
0659 end
0660
0661 Pmax(i) = mpc.gen(genidx(i),PMAX) / Sbase;
0662 Pmin(i) = mpc.gen(genidx(i),PMIN) / Sbase;
0663
0664 if (Pmax(i) - Pmin(i) >= minPgendiff)
0665
0666
0667 constraints = [constraints;
0668 psiL_sdp{k}(i) >= 0;
0669 psiU_sdp{k}(i) >= 0];
0670
0671 if pwl
0672
0673 lambdasum = 0;
0674 for pwiter = 1:length(pwl_sdp{k}{i})
0675
0676 m_pw = (c_pw(pwiter+1) - c_pw(pwiter)) / (x_pw(pwiter+1) - x_pw(pwiter));
0677
0678 if exist('cost','var')
0679 cost = cost - pwl_sdp{k}{i}(pwiter) * (m_pw*x_pw(pwiter) - c_pw(pwiter));
0680 else
0681 cost = -pwl_sdp{k}{i}(pwiter) * (m_pw*x_pw(pwiter) - c_pw(pwiter));
0682 end
0683
0684 lambdasum = lambdasum + pwl_sdp{k}{i}(pwiter) * m_pw;
0685
0686 end
0687
0688 constraints = [constraints;
0689 pwl_sdp{k}{i}(:) >= 0];
0690
0691 if ~exist('lambdai','var')
0692 lambdai = lambdasum + psiU_sdp{k}(i) - psiL_sdp{k}(i);
0693 else
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703 constraints = [constraints;
0704 lambdasum + psiU_sdp{k}(i) - psiL_sdp{k}(i) - binding_lagrange <= lambdai <= lambdasum + psiU_sdp{k}(i) - psiL_sdp{k}(i) + binding_lagrange];
0705 end
0706
0707 constraints = [constraints;
0708 sum(pwl_sdp{k}{i}) == 1];
0709
0710 else
0711 if c2(genidx(i)) > 0
0712
0713
0714
0715 R{k}{i} = sdpvar(2,1);
0716
0717 constraints = [constraints;
0718 [1 R{k}{i}(1);
0719 R{k}{i}(1) R{k}{i}(2)] >= 0];
0720
0721 if ~exist('lambdai','var')
0722 lambdai = c1(genidx(i)) + 2*sqrt(c2(genidx(i)))*R{k}{i}(1) + psiU_sdp{k}(i) - psiL_sdp{k}(i);
0723 else
0724
0725
0726 constraints = [constraints;
0727 c1(genidx(i)) + 2*sqrt(c2(genidx(i)))*R{k}{i}(1) + psiU_sdp{k}(i) - psiL_sdp{k}(i) == lambdai];
0728 end
0729
0730 if exist('cost','var')
0731 cost = cost - R{k}{i}(2);
0732 else
0733 cost = -R{k}{i}(2);
0734 end
0735 else
0736
0737 R{k}{i} = zeros(2,1);
0738
0739 if ~exist('lambdai','var')
0740 lambdai = c1(genidx(i)) + psiU_sdp{k}(i) - psiL_sdp{k}(i);
0741 else
0742
0743
0744 constraints = [constraints;
0745 c1(genidx(i)) + psiU_sdp{k}(i) - psiL_sdp{k}(i) == lambdai];
0746 end
0747
0748 end
0749 end
0750
0751 if exist('cost','var')
0752 cost = cost ...
0753 + psiL_sdp{k}(i)*Pmin(i) - psiU_sdp{k}(i)*Pmax(i);
0754 else
0755 cost = psiL_sdp{k}(i)*Pmin(i) - psiU_sdp{k}(i)*Pmax(i);
0756 end
0757
0758 else
0759
0760
0761
0762
0763 Pavg = (Pmin(i)+Pmax(i))/2;
0764 Pd(k) = Pd(k) - Pavg;
0765
0766
0767 if pwl
0768
0769 pwidx = find(Pavg <= x_pw,1);
0770 if isempty(pwidx)
0771 pwidx = length(x_pw);
0772 end
0773 pwidx = pwidx - 1;
0774
0775 m_pw = (c_pw(pwidx+1) - c_pw(pwidx)) / (x_pw(pwidx+1) - x_pw(pwidx));
0776 fix_gen_cost = m_pw * (Pavg - x_pw(pwidx)) + c_pw(pwidx);
0777
0778 if exist('cost','var')
0779 cost = cost + fix_gen_cost;
0780 else
0781 cost = fix_gen_cost;
0782 end
0783 else
0784 c0(genidx(i)) = c2(genidx(i))*Pavg^2 + c1(genidx(i))*Pavg + c0(genidx(i));
0785 end
0786
0787
0788 end
0789
0790 if ~pwl
0791 if exist('cost','var')
0792 cost = cost + c0(genidx(i));
0793 else
0794 cost = c0(genidx(i));
0795 end
0796 end
0797 end
0798
0799
0800
0801
0802 if exist('lambdai','var')
0803 lambda_aggregate = lambdai;
0804 else
0805 lambda_aggregate = lambdaeq_sdp(lambda_eq(k));
0806 end
0807
0808 lambdaeq_sdp(k) = -lambda_aggregate;
0809
0810 cost = cost + lambda_aggregate*Pd(k);
0811
0812 A = addToA(Yk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, lambda_aggregate, maxclique);
0813
0814
0815
0816
0817
0818
0819
0820 clear gamma_aggregate
0821 if (Qmax(k) - Qmin(k) >= minQgendiff)
0822
0823 if Qmin(k) >= -maxgenlimit
0824
0825 gamma_aggregate = -gammaL_sdp(gamma_ineq(k));
0826
0827 if exist('cost','var')
0828 cost = cost ...
0829 + gammaL_sdp(gamma_ineq(k))*(Qmin(k) - Qd(k));
0830 else
0831 cost = gammaL_sdp(gamma_ineq(k))*(Qmin(k) - Qd(k));
0832 end
0833
0834
0835 constraints = [constraints;
0836 gammaL_sdp(gamma_ineq(k)) >= 0];
0837 end
0838
0839 if Qmax(k) <= maxgenlimit
0840
0841 if exist('gamma_aggregate','var')
0842 gamma_aggregate = gamma_aggregate + gammaU_sdp(gamma_ineq(k));
0843 else
0844 gamma_aggregate = gammaU_sdp(gamma_ineq(k));
0845 end
0846
0847 if exist('cost','var')
0848 cost = cost ...
0849 - gammaU_sdp(gamma_ineq(k))*(Qmax(k) - Qd(k));
0850 else
0851 cost = -gammaU_sdp(gamma_ineq(k))*(Qmax(k) - Qd(k));
0852 end
0853
0854
0855 constraints = [constraints;
0856 gammaU_sdp(gamma_ineq(k)) >= 0];
0857 end
0858
0859 else
0860 gamma_aggregate = -gammaeq_sdp(gamma_eq(k));
0861
0862 if exist('cost','var')
0863 cost = cost ...
0864 + gammaeq_sdp(gamma_eq(k))*((Qmin(k)+Qmax(k))/2 - Qd(k));
0865 else
0866 cost = gammaeq_sdp(gamma_eq(k))*((Qmin(k)+Qmax(k))/2 - Qd(k));
0867 end
0868
0869
0870 end
0871
0872
0873 if (Qmax(k) <= maxgenlimit || Qmin(k) >= -maxgenlimit)
0874 A = addToA(Yk_(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, gamma_aggregate, maxclique);
0875 end
0876
0877
0878 else
0879 error('sdpopf_solver: Invalid bus type');
0880 end
0881
0882
0883
0884
0885 mu_aggregate = -muL_sdp(mu_ineq(k)) + muU_sdp(mu_ineq(k));
0886 A = addToA(Mk(k), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, mu_aggregate, maxclique);
0887
0888
0889 constraints = [constraints;
0890 muL_sdp(mu_ineq(k)) >= 0;
0891 muU_sdp(mu_ineq(k)) >= 0];
0892
0893
0894 cost = cost + muL_sdp(mu_ineq(k))*Vmin(k)^2 - muU_sdp(mu_ineq(k))*Vmax(k)^2;
0895
0896
0897 if verbose >= 2 && mod(k,ndisplay) == 0
0898 fprintf('SDP creation: Bus %i of %i\n',k,nbus);
0899 end
0900
0901 end
0902
0903
0904
0905
0906
0907 nlconstraint = 0;
0908
0909 for i=1:nbranch
0910 if Smax(i) ~= 0 && Smax(i) < maxlinelimit
0911 if upper(mpopt.opf.flow_lim(1)) == 'S'
0912
0913 nlconstraint = nlconstraint + 1;
0914 cost = cost ...
0915 - (Smax(i)^2*Hsdp{nlconstraint}(1,1) + Hsdp{nlconstraint}(2,2) + Hsdp{nlconstraint}(3,3));
0916
0917 constraints = [constraints; Hsdp{nlconstraint} >= 0];
0918
0919 A = addToA(Ylineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,2), maxclique);
0920 A = addToA(Y_lineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,3), maxclique);
0921
0922 nlconstraint = nlconstraint + 1;
0923 cost = cost ...
0924 - (Smax(i)^2*Hsdp{nlconstraint}(1,1) + Hsdp{nlconstraint}(2,2) + Hsdp{nlconstraint}(3,3));
0925
0926 constraints = [constraints; Hsdp{nlconstraint} >= 0];
0927
0928 A = addToA(Ylinetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,2), maxclique);
0929 A = addToA(Y_linetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, 2*Hsdp{nlconstraint}(1,3), maxclique);
0930
0931 elseif upper(mpopt.opf.flow_lim(1)) == 'P'
0932
0933 nlconstraint = nlconstraint + 1;
0934 cost = cost ...
0935 - Smax(i)*Hsdp(nlconstraint);
0936
0937 A = addToA(Ylineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, Hsdp(nlconstraint), maxclique);
0938
0939 nlconstraint = nlconstraint + 1;
0940 cost = cost ...
0941 - Smax(i)*Hsdp(nlconstraint);
0942
0943 A = addToA(Ylinetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, A, Hsdp(nlconstraint), maxclique);
0944
0945 else
0946 error('sdpopf_solver: Invalid line constraint option');
0947 end
0948 end
0949
0950 if verbose >= 2 && mod(i,ndisplay) == 0
0951 fprintf('SDP creation: Branch %i of %i\n',i,nbranch);
0952 end
0953
0954 end
0955
0956 if upper(mpopt.opf.flow_lim(1)) == 'P'
0957 constraints = [constraints; Hsdp >= 0];
0958 end
0959
0960
0961
0962
0963 Aconstraints = zeros(nmaxclique,1);
0964 for i=1:nmaxclique
0965
0966
0967 constraints = [constraints; 1*A{i} >= 0];
0968 Aconstraints(i) = length(constraints);
0969 end
0970
0971 tsdpform = toc;
0972
0973
0974
0975 if recover_voltage == 2 || recover_voltage == 3 || recover_voltage == 4
0976 sdpopts = sdpsettings(sdpopts,'saveduals',1);
0977 end
0978
0979
0980 S = warning;
0981
0982
0983
0984
0985 sdpinfo = solvesdp(constraints, -cost, sdpopts);
0986
0987 if sdpinfo.problem == 2 || sdpinfo.problem == -3
0988 error(yalmiperror(sdpinfo.problem));
0989 end
0990 warning(S);
0991 tsolve = sdpinfo.solvertime;
0992
0993 if verbose >= 2
0994 fprintf('Solver exit message: %s\n',sdpinfo.info);
0995 end
0996
0997 objval = double(cost);
0998
0999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009 tic;
1010
1011
1012 if recover_voltage == 1 || recover_voltage == 3 || recover_voltage == 4
1013 dual_mineigratio = inf;
1014 dual_eval = nan(length(A),1);
1015 for i=1:length(A)
1016 [evc, evl] = eig(double(A{i}));
1017
1018
1019 evl = diag(evl);
1020 dual_u{i} = evc(:,abs(evl) == min(abs(evl)));
1021
1022
1023 dual_eigA_ratio = abs(evl(3) / evl(1));
1024 dual_eval(i) = 1/dual_eigA_ratio;
1025
1026 if dual_eigA_ratio < dual_mineigratio
1027 dual_mineigratio = dual_eigA_ratio;
1028 end
1029 end
1030 end
1031
1032
1033 if recover_voltage == 2 || recover_voltage == 3 || recover_voltage == 4
1034 primal_mineigratio = inf;
1035 primal_eval = nan(length(A),1);
1036 for i=1:length(A)
1037 [evc, evl] = eig(double(dual(constraints(Aconstraints(i)))));
1038
1039
1040 evl = diag(evl);
1041 primal_u{i} = evc(:,abs(evl) == max(abs(evl)));
1042
1043
1044 primal_eigA_ratio = abs(evl(end) / evl(end-2));
1045 primal_eval(i) = 1/primal_eigA_ratio;
1046
1047 if primal_eigA_ratio < primal_mineigratio
1048 primal_mineigratio = primal_eigA_ratio;
1049 end
1050 end
1051 end
1052
1053 if recover_voltage == 3 || recover_voltage == 4
1054 recover_voltage_loop = [1 2];
1055 else
1056 recover_voltage_loop = recover_voltage;
1057 end
1058
1059 for r = 1:length(recover_voltage_loop)
1060
1061 if recover_voltage_loop(r) == 1
1062 if verbose >= 2
1063 fprintf('Recovering dual solution.\n');
1064 end
1065 u = dual_u;
1066 eval = dual_eval;
1067 mineigratio = dual_mineigratio;
1068 elseif recover_voltage_loop(r) == 2
1069 if verbose >= 2
1070 fprintf('Recovering primal solution.\n');
1071 end
1072 u = primal_u;
1073 eval = primal_eval;
1074 mineigratio = primal_mineigratio;
1075 else
1076 error('sdpopf_solver: Invalid recover_voltage');
1077 end
1078
1079
1080
1081
1082
1083 Aint = [];
1084 Aslack = [];
1085 slackbus_idx = find(mpc.bus(:,2) == 3);
1086 for i=1:nmaxclique
1087 maxcliquei = maxclique{i};
1088 for k=i+1:nmaxclique
1089
1090 temp = maxcliquei(ismembc(maxcliquei, maxclique{k}));
1091
1092
1093
1094
1095
1096 Aint = [Aint; i*ones(length(temp),1) k*ones(length(temp),1) temp(:)];
1097 end
1098
1099 if any(maxcliquei == slackbus_idx)
1100 Aslack = i;
1101 end
1102 end
1103 if nmaxclique > 1
1104 L = sparse(2*length(Aint(:,1))+1,2*nmaxclique);
1105 else
1106 L = sparse(1,2);
1107 end
1108 Lrow = 0;
1109 for i=1:size(Aint,1)
1110
1111 Lrow = Lrow + 1;
1112 L(Lrow,Aint(i,1)) = u{Aint(i,1)}(maxclique{Aint(i,1)}(:) == Aint(i,3));
1113 L(Lrow,Aint(i,1)+nmaxclique) = -u{Aint(i,1)}(find(maxclique{Aint(i,1)}(:) == Aint(i,3))+length(maxclique{Aint(i,1)}(:)));
1114 L(Lrow,Aint(i,2)) = -u{Aint(i,2)}(maxclique{Aint(i,2)}(:) == Aint(i,3));
1115 L(Lrow,Aint(i,2)+nmaxclique) = u{Aint(i,2)}(find(maxclique{Aint(i,2)}(:) == Aint(i,3))+length(maxclique{Aint(i,2)}(:)));
1116
1117
1118 Lrow = Lrow + 1;
1119 L(Lrow,Aint(i,1)) = u{Aint(i,1)}(find(maxclique{Aint(i,1)}(:) == Aint(i,3))+length(maxclique{Aint(i,1)}(:)));
1120 L(Lrow,Aint(i,1)+nmaxclique) = u{Aint(i,1)}(maxclique{Aint(i,1)}(:) == Aint(i,3));
1121 L(Lrow,Aint(i,2)) = -u{Aint(i,2)}(find(maxclique{Aint(i,2)}(:) == Aint(i,3))+length(maxclique{Aint(i,2)}(:)));
1122 L(Lrow,Aint(i,2)+nmaxclique) = -u{Aint(i,2)}(maxclique{Aint(i,2)}(:) == Aint(i,3));
1123 end
1124
1125
1126 Lrow = Lrow + 1;
1127 L(Lrow,Aslack) = u{Aslack}(find(maxclique{Aslack}(:) == slackbus_idx)+length(maxclique{Aslack}(:)));
1128 L(Lrow,Aslack+nmaxclique) = u{Aslack}(maxclique{Aslack}(:) == slackbus_idx);
1129
1130
1131 if nmaxclique == 1
1132 [LLevec,LLeval] = eig(full(L).'*full(L));
1133 LLeval = LLeval(1);
1134 LLevec = LLevec(:,1);
1135 else
1136
1137
1138
1139
1140
1141 try
1142 warning off
1143 [LLevec,LLeval] = eigs(L.'*L,1,'SM');
1144 warning(S);
1145 catch
1146 [LLevec,LLeval] = eig(full(L).'*full(L));
1147 LLeval = LLeval(1);
1148 LLevec = LLevec(:,1);
1149 warning(S);
1150 end
1151 end
1152
1153
1154
1155
1156
1157
1158 U = sparse(2*nbus,nmaxclique);
1159 Ueval = sparse(2*nbus,nmaxclique);
1160 for i=1:nmaxclique
1161 newentries = [maxclique{i}(:); maxclique{i}(:)+nbus];
1162 newvals_phasor = (u{i}(1:length(maxclique{i}))+1i*u{i}(length(maxclique{i})+1:2*length(maxclique{i})))*(LLevec(i)+1i*LLevec(i+nmaxclique));
1163 newvals = [real(newvals_phasor(:)); imag(newvals_phasor(:))];
1164
1165 Ueval(newentries,i) = 1/eval(i);
1166 U(newentries,i) = newvals;
1167
1168 end
1169 Ueval_sum = sum(Ueval,2);
1170 for i=1:2*nbus
1171 Ueval(i,:) = Ueval(i,:) / Ueval_sum(i);
1172 end
1173 U = full(sum(U .* Ueval,2));
1174
1175
1176 muL_val = abs(double(muL_sdp(mu_ineq)));
1177 muU_val = abs(double(muU_sdp(mu_ineq)));
1178
1179 [maxmuU,binding_voltage_busU] = max(muU_val);
1180 [maxmuL,binding_voltage_busL] = max(muL_val);
1181
1182 if (maxmuU < binding_lagrange) && (maxmuL < binding_lagrange)
1183 if verbose > 1
1184 fprintf('No binding voltage magnitude. Attempting solution recovery from line-flow constraints.\n');
1185 end
1186
1187
1188
1189
1190
1191
1192 largest_H_norm = 0;
1193 largest_H_norm_idx = 0;
1194 for i=1:length(Hsdp)
1195 if norm(double(Hsdp{i})) > largest_H_norm
1196 largest_H_norm = norm(double(Hsdp{i}));
1197 largest_H_norm_idx = i;
1198 end
1199 end
1200
1201 if largest_H_norm < binding_lagrange
1202 if verbose > 1
1203 error('sdpopf_solver: No binding voltage magnitude or line-flow limits, and no PQ buses. Error recovering solution.');
1204 end
1205 else
1206
1207 if Hlookup(largest_H_norm_idx,2) == 1
1208 [Yline_row, Yline_col, Yline_val] = find(Ylineft(Hlookup(largest_H_norm_idx,1)));
1209 [Y_line_row, Y_line_col, Y_line_val] = find(Y_lineft(Hlookup(largest_H_norm_idx,1)));
1210 else
1211 [Yline_row, Yline_col, Yline_val] = find(Ylinetf(Hlookup(largest_H_norm_idx,1)));
1212 [Y_line_row, Y_line_col, Y_line_val] = find(Y_linetf(Hlookup(largest_H_norm_idx,1)));
1213 end
1214 trYlineUUT = 0;
1215 for i=1:length(Yline_val)
1216 trYlineUUT = trYlineUUT + Yline_val(i) * U(Yline_row(i))*U(Yline_col(i));
1217 end
1218
1219 if upper(mpopt.opf.flow_lim(1)) == 'S'
1220 trY_lineUUT = 0;
1221 for i=1:length(Y_line_val)
1222 trY_lineUUT = trY_lineUUT + Y_line_val(i) * U(Y_line_row(i))*U(Y_line_col(i));
1223 end
1224 chi = (Smax(Hlookup(largest_H_norm_idx,1))^2 / ( trYlineUUT^2 + trY_lineUUT^2 ) )^(0.25);
1225 elseif upper(mpopt.opf.flow_lim(1)) == 'P'
1226
1227 chi = sqrt(Smax(Hlookup(largest_H_norm_idx,1)) / trYlineUUT);
1228 end
1229
1230 U = chi*U;
1231 end
1232
1233 else
1234 if maxmuU > maxmuL
1235 binding_voltage_bus = binding_voltage_busU;
1236 binding_voltage_mag = Vmax(binding_voltage_bus);
1237 else
1238 binding_voltage_bus = binding_voltage_busL;
1239 binding_voltage_mag = Vmin(binding_voltage_bus);
1240 end
1241 chi = binding_voltage_mag^2 / (U(binding_voltage_bus)^2 + U(binding_voltage_bus+nbus)^2);
1242 U = sqrt(chi)*U;
1243 end
1244
1245 Vopt = (U(1:nbus) + 1i*U(nbus+1:2*nbus));
1246
1247 if real(Vopt(slackbus_idx)) < 0
1248 Vopt = -Vopt;
1249 end
1250
1251 if recover_voltage == 3 || recover_voltage == 4
1252 Sopt = Vopt .* conj(Y*Vopt);
1253 PQerr = norm([real(Sopt(bustype(:) == PQ))*100+mpc.bus(bustype(:) == PQ,PD); imag(Sopt(bustype(:) == PQ))*100+mpc.bus(bustype(:) == PQ,QD)], 2);
1254 if recover_voltage_loop(r) == 1
1255 Vopt_dual = Vopt;
1256 dual_PQerr = PQerr;
1257 elseif recover_voltage_loop(r) == 2
1258 Vopt_primal = Vopt;
1259 primal_PQerr = PQerr;
1260 end
1261 end
1262 end
1263
1264 if recover_voltage == 3
1265 if dual_PQerr > primal_PQerr
1266 fprintf('dual_PQerr: %g MW. primal_PQerr: %g MW. Using primal solution.\n',dual_PQerr,primal_PQerr)
1267 Vopt = Vopt_primal;
1268 else
1269 fprintf('dual_PQerr: %g MW. primal_PQerr: %g MW. Using dual solution.\n',dual_PQerr,primal_PQerr)
1270 Vopt = Vopt_dual;
1271 end
1272 elseif recover_voltage == 4
1273 if primal_mineigratio > dual_mineigratio
1274 if verbose >= 2
1275 fprintf('dual_mineigratio: %g. primal_mineigratio: %g. Using primal solution.\n',dual_mineigratio,primal_mineigratio)
1276 end
1277 Vopt = Vopt_primal;
1278 else
1279 if verbose >= 2
1280 fprintf('dual_mineigratio: %g. primal_mineigratio: %g. Using dual solution.\n',dual_mineigratio,primal_mineigratio)
1281 end
1282 Vopt = Vopt_dual;
1283 end
1284 end
1285
1286
1287 for i=1:nbus
1288 mpc.bus(i,VM) = abs(Vopt(i));
1289 mpc.bus(i,VA) = angle(Vopt(i))*180/pi;
1290 mpc.gen(mpc.gen(:,GEN_BUS) == i,VG) = abs(Vopt(i));
1291 end
1292
1293
1294
1295
1296 if recover_injections == 2
1297
1298
1299 for i=1:length(A)
1300 W{i} = dual(constraints(Aconstraints(i)));
1301 end
1302
1303 Pinj = zeros(nbus,1);
1304 Qinj = zeros(nbus,1);
1305 for i=1:nbus
1306
1307
1308 Pinj(i) = recoverFromW(Yk(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique);
1309
1310
1311 Qinj(i) = recoverFromW(Yk_(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique);
1312
1313
1314
1315
1316 end
1317
1318 Pflowft = zeros(nbranch,1);
1319 Pflowtf = zeros(nbranch,1);
1320 Qflowft = zeros(nbranch,1);
1321 Qflowtf = zeros(nbranch,1);
1322 for i=1:nbranch
1323 Pflowft(i) = recoverFromW(Ylineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1324 Pflowtf(i) = recoverFromW(Ylinetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1325 Qflowft(i) = recoverFromW(Y_lineft(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1326 Qflowtf(i) = recoverFromW(Y_linetf(i), Wref_dd, Wref_qq, Wref_dq, matidx_dd, matidx_qq, matidx_dq, W, maxclique) * Sbase;
1327 end
1328
1329
1330 elseif recover_injections == 1
1331
1332 Sopt = Vopt .* conj(Y*Vopt);
1333 Pinj = real(Sopt);
1334 Qinj = imag(Sopt);
1335
1336 Vf = Vopt(mpc.branch(:,F_BUS));
1337 Vt = Vopt(mpc.branch(:,T_BUS));
1338
1339 Sft = Vf.*conj(Yf*Vopt);
1340 Stf = Vt.*conj(Yt*Vopt);
1341
1342 Pflowft = real(Sft) * Sbase;
1343 Pflowtf = real(Stf) * Sbase;
1344 Qflowft = imag(Sft) * Sbase;
1345 Qflowtf = imag(Stf) * Sbase;
1346 end
1347
1348
1349 busorder = mpc.bus(:,end);
1350 mpc.bus = mpc.bus(:,1:end-1);
1351
1352
1353 Qtol = 5e-2;
1354 mpc.gen(:,PG) = 0;
1355 mpc.gen(:,QG) = 0;
1356 for i=1:nbus
1357
1358 if bustype(i) == PV || bustype(i) == REF
1359
1360 genidx = find(mpc.gen(:,GEN_BUS) == i);
1361
1362
1363
1364
1365
1366
1367
1368 gen_not_limited_pw = [];
1369 gen_not_limited_q = [];
1370 Pmax = mpc.gen(genidx,PMAX);
1371 Pmin = mpc.gen(genidx,PMIN);
1372 remove_from_list = [];
1373 for k=1:length(genidx)
1374
1375 if (Pmax(k) - Pmin(k) < minPgendiff * Sbase)
1376 mpc.gen(genidx(k),PG) = (Pmin(k)+Pmax(k))/2;
1377 remove_from_list = [remove_from_list; k];
1378 else
1379 if mpc.gencost(genidx(k),MODEL) == PW_LINEAR
1380 gen_not_limited_pw = [gen_not_limited_pw; genidx(k)];
1381 mpc.gen(genidx(k),PG) = mpc.gen(genidx(k),PMIN);
1382 else
1383 gen_not_limited_q = [gen_not_limited_q; genidx(k)];
1384 mpc.gen(genidx(k),PG) = 0;
1385 end
1386 end
1387 end
1388 Pmax(remove_from_list) = [];
1389 Pmin(remove_from_list) = [];
1390
1391
1392
1393 Pgen_remain = Pinj(i) + mpc.bus(i,PD) / Sbase - sum(mpc.gen(genidx,PG)) / Sbase;
1394
1395
1396 gen_not_limited = [gen_not_limited_q; gen_not_limited_pw];
1397 if length(gen_not_limited) == 1
1398
1399
1400 mpc.gen(gen_not_limited,PG) = mpc.gen(gen_not_limited,PG) + Pgen_remain * Sbase;
1401 elseif length(genidx) == 1
1402
1403
1404 mpc.gen(genidx,PG) = mpc.gen(genidx,PG) + Pgen_remain * Sbase;
1405 else
1406
1407
1408
1409
1410
1411
1412
1413 if ~isempty(gen_not_limited_pw)
1414 while Pgen_remain > 1e-5
1415
1416
1417 lc = inf;
1418 for geniter = 1:length(gen_not_limited_pw)
1419 nseg = mpc.gencost(gen_not_limited_pw(geniter),NCOST)-1;
1420
1421
1422 x_pw = mpc.gencost(gen_not_limited_pw(geniter),-1+COST+(1:2:2*nseg+1)) / Sbase;
1423 c_pw = mpc.gencost(gen_not_limited_pw(geniter),-1+COST+(2:2:2*nseg+2));
1424
1425 for pwidx = 1:nseg
1426 m_pw = (c_pw(pwidx+1) - c_pw(pwidx)) / (x_pw(pwidx+1) - x_pw(pwidx));
1427 if x_pw(pwidx) > mpc.gen(gen_not_limited_pw(geniter),PG)/Sbase
1428 x_min = x_pw(pwidx-1);
1429 x_max = x_pw(pwidx);
1430 break;
1431 end
1432 end
1433
1434 if m_pw < lc
1435 lc = m_pw;
1436 lc_x_min = x_min;
1437 lc_x_max = x_max;
1438 lc_geniter = geniter;
1439 end
1440
1441 end
1442 if mpc.gen(gen_not_limited_pw(lc_geniter), PG) + min(Pgen_remain, lc_x_max-lc_x_min)*Sbase < Pmax(lc_geniter)
1443
1444
1445 mpc.gen(gen_not_limited_pw(lc_geniter), PG) = mpc.gen(gen_not_limited_pw(lc_geniter), PG) + min(Pgen_remain, lc_x_max-lc_x_min)*Sbase;
1446 Pgen_remain = Pgen_remain - min(Pgen_remain, lc_x_max-lc_x_min);
1447 else
1448
1449
1450
1451 Pgen_remain = Pgen_remain - (Pmax(lc_geniter) - mpc.gen(gen_not_limited_pw(lc_geniter), PG));
1452 mpc.gen(gen_not_limited_pw(lc_geniter), PG) = Pmax(lc_geniter);
1453 gen_not_limited_pw(lc_geniter) = [];
1454 Pmax(lc_geniter) = [];
1455 Pmin(lc_geniter) = [];
1456 end
1457 end
1458 end
1459
1460
1461
1462
1463
1464 if ~isempty(gen_not_limited_q)
1465
1466 Pg_q = sdpvar(length(gen_not_limited_q),1);
1467 cost_q = Pg_q.' * diag(c2(gen_not_limited_q)) * Pg_q + c1(gen_not_limited_q).' * Pg_q;
1468
1469
1470
1471
1472
1473 if isnumeric(cost_q)
1474 cost_q = sum(Pg_q);
1475 end
1476
1477 constraints_q = [Pmin/Sbase <= Pg_q <= Pmax/Sbase;
1478 sum(Pg_q) == Pgen_remain];
1479 sol_q = solvesdp(constraints_q,cost_q, sdpsettings(sdpopts,'Verbose',0,'saveduals',0));
1480
1481
1482
1483
1484
1485
1486
1487
1488 constraint_q_err = checkset(constraints_q);
1489 if sol_q.problem == 0 || max(abs(constraint_q_err)) < 1e-3
1490 mpc.gen(gen_not_limited_q, PG) = double(Pg_q) * Sbase;
1491 elseif sol_q.problem == 4 || sol_q.problem == 1
1492
1493
1494
1495 mpc.gen(gen_not_limited_q, PG) = double(Pg_q) * Sbase;
1496 warning('sdpopf_solver: Numeric inconsistency assigning active power injection to generators at bus %i.',busorder(i));
1497 else
1498
1499
1500
1501 mpc.gen(gen_not_limited_q, PG) = (Pgen_remain / length(gen_not_limited_q))*Sbase;
1502 warning('sdpopf_solver: Error assigning active power injection to generators at bus %i. Assiging equally among all generators at this bus.',busorder(i));
1503 end
1504
1505 end
1506
1507 end
1508
1509
1510
1511
1512 mpc.gen(genidx,QG) = mpc.gen(genidx,QMIN);
1513 Qremaining = (Qinj(i) + Qd(i))*Sbase - sum(mpc.gen(genidx,QG));
1514 genidx_idx = 1;
1515 while abs(Qremaining) > Qtol && genidx_idx <= length(genidx)
1516 if mpc.gen(genidx(genidx_idx),QG) + Qremaining <= mpc.gen(genidx(genidx_idx),QMAX)
1517 mpc.gen(genidx(genidx_idx),QG) = mpc.gen(genidx(genidx_idx),QG) + Qremaining;
1518 Qremaining = 0;
1519 else
1520 Qremaining = Qremaining - (mpc.gen(genidx(genidx_idx),QMAX) - mpc.gen(genidx(genidx_idx),QG));
1521 mpc.gen(genidx(genidx_idx),QG) = mpc.gen(genidx(genidx_idx),QMAX);
1522 end
1523 genidx_idx = genidx_idx + 1;
1524 end
1525 if abs(Qremaining) > Qtol
1526
1527
1528
1529 mpc.gen(genidx,QG) = mpc.gen(genidx,QG) + Qremaining / length(genidx);
1530 warning('sdpopf_solver: Inconsistency in assigning reactive power to generators at bus index %i.',i);
1531 end
1532 end
1533 end
1534
1535
1536 mpc.branch(:,PF) = Pflowft;
1537 mpc.branch(:,PT) = Pflowtf;
1538 mpc.branch(:,QF) = Qflowft;
1539 mpc.branch(:,QT) = Qflowtf;
1540
1541
1542
1543 mpc.bus(:,LAM_P) = -double(lambdaeq_sdp) / Sbase;
1544 for i=1:nbus
1545 if ~isnan(gamma_eq(i))
1546 mpc.bus(i,LAM_Q) = -double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1547 elseif Qmin(i) >= -maxgenlimit && Qmax(i) <= maxgenlimit
1548 mpc.bus(i,LAM_Q) = (double(gammaU_sdp(gamma_ineq(i))) - double(gammaL_sdp(gamma_ineq(i)))) / Sbase;
1549 elseif Qmax(i) <= maxgenlimit
1550 mpc.bus(i,LAM_Q) = double(gammaU_sdp(gamma_ineq(i))) / Sbase;
1551 elseif Qmin(i) >= -maxgenlimit
1552 mpc.bus(i,LAM_Q) = double(gammaL_sdp(gamma_ineq(i))) / Sbase;
1553 else
1554 mpc.bus(i,LAM_Q) = 0;
1555 end
1556 end
1557
1558
1559
1560
1561 mpc.bus(:,MU_VMAX) = 2*abs(Vopt) .* double(muU_sdp);
1562 mpc.bus(:,MU_VMIN) = 2*abs(Vopt) .* double(muL_sdp);
1563
1564
1565 for i=1:nbus
1566 genidx = find(mpc.gen(:,GEN_BUS) == i);
1567 for k=1:length(genidx)
1568 if ~isnan(double(psiU_sdp{i}(k))) && ~isnan(double(psiU_sdp{i}(k)))
1569 mpc.gen(genidx(k),MU_PMAX) = double(psiU_sdp{i}(k)) / Sbase;
1570 mpc.gen(genidx(k),MU_PMIN) = double(psiL_sdp{i}(k)) / Sbase;
1571 else
1572
1573
1574
1575 Pavg = (mpc.gen(genidx(k),PMAX)+mpc.gen(genidx(k),PMIN))/(2*Sbase);
1576 c0k = mpc.gencost(genidx(k),COST+2);
1577 alphak = c2(genidx(k))*Pavg^2 + c1(genidx(k))*Pavg + c0k;
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590 a = 1;
1591 b = -2*sqrt(c2(genidx(k)))*Pavg;
1592 c = -c1(genidx(k))*Pavg+alphak-c0k;
1593 R12 = (-b+sqrt(b^2-4*a*c)) / (2*a);
1594 if abs(imag(R12)) > 1e-3
1595 warning('sdpopf_solver: R12 has a nonzero imaginary component. Lagrange multipliers for generator active power limits may be incorrect.');
1596 end
1597 R12 = real(R12);
1598
1599
1600
1601 psi_tight = -double(lambdaeq_sdp(lambda_eq(i))) / Sbase - c1(genidx(k)) / Sbase - 2*sqrt(c2(genidx(k))/Sbase^2)*R12;
1602 if psi_tight < 0
1603 mpc.gen(genidx(k),MU_PMAX) = 0;
1604 mpc.gen(genidx(k),MU_PMIN) = -psi_tight;
1605 else
1606 mpc.gen(genidx(k),MU_PMAX) = psi_tight;
1607 mpc.gen(genidx(k),MU_PMIN) = 0;
1608 end
1609 end
1610 end
1611 if ~isempty(genidx)
1612 if ~isnan(gamma_eq(i))
1613 gamma_tight = double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1614 if gamma_tight > 0
1615 mpc.gen(genidx(k),MU_QMAX) = 0;
1616 mpc.gen(genidx(k),MU_QMIN) = double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1617 else
1618 mpc.gen(genidx(k),MU_QMAX) = -double(gammaeq_sdp(gamma_eq(i))) / Sbase;
1619 mpc.gen(genidx(k),MU_QMIN) = 0;
1620 end
1621 else
1622 if any(Qmin(i) < -maxgenlimit)
1623 mpc.gen(genidx,MU_QMIN) = 0;
1624 else
1625 mpc.gen(genidx,MU_QMIN) = double(gammaL_sdp(gamma_ineq(i))) / Sbase;
1626 end
1627
1628 if any(Qmax(i) > maxgenlimit)
1629 mpc.gen(genidx,MU_QMAX) = 0;
1630 else
1631 mpc.gen(genidx,MU_QMAX) = double(gammaU_sdp(gamma_ineq(i))) / Sbase;
1632 end
1633 end
1634 end
1635 end
1636
1637
1638
1639
1640
1641
1642
1643 if nlconstraint > 0
1644 for i=1:nbranch
1645 Hidx = Hlookup(:,1) == i & Hlookup(:,2) == 1;
1646 if any(Hidx)
1647 if upper(mpopt.opf.flow_lim(1)) == 'S'
1648 Hft = double(Hsdp{Hidx}) ./ Sbase;
1649 mpc.branch(i,MU_SF) = 2*sqrt((Hft(1,2)^2 + Hft(1,3)^2));
1650 Htf = double(Hsdp{Hlookup(:,1) == i & Hlookup(:,2) == 0}) ./ Sbase;
1651 mpc.branch(i,MU_ST) = 2*sqrt((Htf(1,2)^2 + Htf(1,3)^2));
1652 elseif upper(mpopt.opf.flow_lim(1)) == 'P'
1653 mpc.branch(i,MU_SF) = Hsdp(Hidx) / Sbase;
1654 mpc.branch(i,MU_ST) = Hsdp(Hlookup(:,1) == i & Hlookup(:,2) == 0) / Sbase;
1655 end
1656 else
1657 mpc.branch(i,MU_SF) = 0;
1658 mpc.branch(i,MU_ST) = 0;
1659 end
1660 end
1661 else
1662 mpc.branch(:,MU_SF) = 0;
1663 mpc.branch(:,MU_ST) = 0;
1664 end
1665
1666
1667 mpc.branch(:,MU_ANGMIN) = nan;
1668 mpc.branch(:,MU_ANGMAX) = nan;
1669
1670
1671 mpc.f = objval;
1672
1673
1674
1675 mpc.bus(:,BUS_I) = busorder;
1676 mpc.bus = sortrows(mpc.bus,1);
1677
1678 for i=1:ngen
1679 mpc.gen(i,GEN_BUS) = busorder(mpc.gen(i,GEN_BUS));
1680 end
1681
1682 [mpc.gen genidx] = sortrows(mpc.gen,1);
1683 mpc.gencost = mpc.gencost(genidx,:);
1684
1685 for i=1:nbranch
1686 mpc.branch(i,F_BUS) = busorder(mpc.branch(i,F_BUS));
1687 mpc.branch(i,T_BUS) = busorder(mpc.branch(i,T_BUS));
1688 end
1689
1690
1691
1692
1693 results = mpc;
1694 success = mineigratio > mineigratio_tol && LLeval < zeroeval_tol;
1695
1696 raw.xr.A = A;
1697 if recover_injections == 2
1698 raw.xr.W = W;
1699 else
1700 raw.xr.W = [];
1701 end
1702
1703 raw.pimul.lambdaeq_sdp = lambdaeq_sdp;
1704 raw.pimul.gammaeq_sdp = gammaeq_sdp;
1705 raw.pimul.gammaU_sdp = gammaU_sdp;
1706 raw.pimul.gammaL_sdp = gammaL_sdp;
1707 raw.pimul.muU_sdp = muU_sdp;
1708 raw.pimul.muL_sdp = muL_sdp;
1709 raw.pimul.psiU_sdp = psiU_sdp;
1710 raw.pimul.psiL_sdp = psiL_sdp;
1711
1712 raw.info = sdpinfo.problem;
1713
1714 results.zero_eval = LLeval;
1715 results.mineigratio = mineigratio;
1716
1717 results.sdpinfo = sdpinfo;
1718
1719 toutput = toc;
1720 results.et = tsetup + tmaxclique + tsdpform + tsolve + toutput;