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