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