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