0001 function [Pbus,Qbus] = sgvm_calc_injection_delta(mpc, opt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 define_constants;
0022 prob = struct();
0023
0024 if isfield(opt.vm.nodeperm, 'nox')
0025 nox = opt.vm.nodeperm.nox;
0026 else
0027 nox = true;
0028 end
0029
0030 if isfield(opt.vm.nodeperm, 'usedv')
0031 usedv = opt.vm.nodeperm.usedv;
0032 else
0033 usedv = false;
0034 end
0035
0036 if isfield(opt.vm.nodeperm, 'branch_slack')
0037 branch_slack = opt.vm.nodeperm.branch_slack;
0038 else
0039 branch_slack = false;
0040 end
0041
0042 if isfield(opt.vm.nodeperm, 'scale_s')
0043 scale_s = opt.vm.nodeperm.scale_s;
0044 else
0045 scale_s = 1;
0046 end
0047
0048 nb = size(mpc.bus,1);
0049 nl = size(mpc.branch,1);
0050
0051
0052 if usedv
0053
0054 Ybus = makeYbus(mpc);
0055 V = mpc.bus(:,VM).*exp(1i*mpc.bus(:,VA)*pi/180);
0056 [~, dsbus_dvm] = dSbus_dV(Ybus,V);
0057 Jpv = real(dsbus_dvm);
0058 Jqv = imag(dsbus_dvm);
0059 clear dsbus_dvm;
0060
0061
0062 Jpv(abs(Jpv) < 1e-6) = 0;
0063 Jqv(abs(Jqv) < 1e-6) = 0;
0064
0065 [rjpv, cjpv, vjpv] = find(Jpv);
0066 [rjqv, cjqv, vjqv] = find(Jqv);
0067 clear Jpv Jqv;
0068
0069 V0 = abs(V);
0070 Vmin = mpc.bus(:,VMIN);
0071 Vmax = mpc.bus(:,VMAX);
0072 end
0073 Sbus = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
0074 dPmax = max(real(Sbus)) - min(real(Sbus));
0075 dQmax = max(imag(Sbus)) - min(imag(Sbus));
0076 wsv = max(abs(Sbus));
0077
0078
0079 Pl0 = mpc.branch(:,PF) / mpc.baseMVA;
0080 Ql0 = mpc.branch(:,QF) / mpc.baseMVA;
0081 Smax = scale_s * mpc.branch(:,RATE_A) / mpc.baseMVA;
0082 wspq = nb*max(Smax);
0083
0084
0085 Smax(Smax == 0) = 5*max(Smax);
0086
0087
0088 alpha = atan(Ql0./Pl0);
0089 alpha(Pl0 < 0) = pi - alpha(Pl0 < 0 );
0090
0091
0092 lfrac = sort(sqrt(Pl0.^2 + Ql0.^2)./Smax, 'descend');
0093 if (nl < 1000) || (lfrac(1000) < 0.5)
0094 nluse = find(sqrt(Pl0.^2 + Ql0.^2)./Smax > 0.5);
0095 else
0096 nluse = find(sqrt(Pl0.^2 + Ql0.^2)./Smax > 0.75);
0097 end
0098 nl = length(nluse);
0099
0100
0101 tptdf = tic;
0102 H = sgvm_acptdf(mpc, nluse);
0103 [rxpp, cxpp, vxpp] = find(H(1:nl, 1:nb));
0104 [rxqp, cxqp, vxqp] = find(H(1:nl, nb + 1: 2*nb));
0105 [rxpq, cxpq, vxpq] = find(H(nl+1:end, 1:nb));
0106 [rxqq, cxqq, vxqq] = find(H(nl+1:end, nb + 1: 2*nb));
0107
0108
0109
0110
0111 clear H;
0112 if opt.vm.nodeperm.verbose > 1
0113 fprintf('\t AC-PTDF time %0.3f sec\n', toc(tptdf));
0114 end
0115
0116
0117 Smax = Smax(nluse);
0118 Pl0 = Pl0(nluse);
0119 Ql0 = Ql0(nluse);
0120 alpha = alpha(nluse);
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132 tsetup = tic;
0133 vars = struct('xp', struct('first', 1, 'last', nb),...
0134 'xq', struct('first',nb+1, 'last',2*nb),...
0135 'absx', struct('first',2*nb+1, 'last',4*nb));
0136 ptr = 4*nb;
0137 if branch_slack
0138 vars.sp = struct('first', ptr+1, 'last', ptr+nl);
0139 ptr = ptr + nl;
0140 vars.sq = struct('first', ptr+1, 'last', ptr+nl);
0141 ptr = ptr + nl;
0142 end
0143 if usedv
0144 vars.dv = struct('first',ptr+1, 'last',ptr + nb);
0145 ptr = ptr + nb;
0146 vars.sv = struct('first', ptr+1, 'last', ptr+1);
0147 ptr = ptr + 1;
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165 end
0166 total_vars = ptr;
0167 prob.A = [];
0168 prob.l = [];
0169 prob.u = [];
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182 if usedv
0183 if ~nox
0184 A = sparse([1:nb, rjpv.' ], ...
0185 [vars.xp.first:vars.xp.last, vars.dv.first - 1 + cjpv.'],...
0186 [ones(1,nb), vjpv.'], nb, total_vars);
0187 lb = zeros(nb, 1);
0188 ub = zeros(nb, 1);
0189 prob = update_Albub(prob, A, lb, ub);
0190 end
0191
0192 A = sparse([1:nb, rjqv.' ], ...
0193 [vars.xq.first:vars.xq.last, vars.dv.first - 1 + cjqv.'],...
0194 [ones(1,nb), vjqv.'], nb, total_vars);
0195 lb = zeros(nb, 1);
0196 ub = zeros(nb, 1);
0197 prob = update_Albub(prob, A, lb, ub);
0198 end
0199
0200
0201
0202
0203
0204
0205
0206
0207 if usedv
0208 A = sparse([1:nb, 1:nb],...
0209 [vars.dv.first:vars.dv.last, vars.sv.first*ones(1,nb)],...
0210 [ones(1,nb), ones(1,nb)], nb, total_vars);
0211
0212
0213
0214 lb = Vmin - V0;
0215 lb(abs(lb) < 1e-6) = 0;
0216 ub = Inf(nb,1);
0217 prob = update_Albub(prob, A, lb, ub);
0218
0219 A = sparse([1:nb, 1:nb],...
0220 [vars.dv.first:vars.dv.last, vars.sv.first*ones(1,nb)],...
0221 [ones(1,nb), -ones(1,nb)], nb, total_vars);
0222
0223
0224
0225
0226 lb = -Inf(nb,1);
0227 ub = Vmax - V0;
0228 ub(abs(ub) < 1e-6) = 0;
0229 prob = update_Albub(prob, A, lb, ub);
0230 end
0231
0232
0233
0234
0235
0236
0237
0238
0239 r = rxpp.';
0240 c = vars.xp.first - 1 + cxpp.';
0241 v = vxpp.';
0242 if ~nox
0243 r = [r, rxqp.'];
0244 c = [c, vars.xq.first - 1 + cxqp.'];
0245 v = [v, vxqp.'];
0246 end
0247 if branch_slack
0248 r = [r, 1:nl];
0249 c = [c, vars.sp.first:vars.sp.last];
0250 v = [v, ones(1,nl)];
0251 end
0252 A = sparse(r, c, v, nl, total_vars);
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264 lb = -abs(Smax.*cos(alpha)) - Pl0;
0265 lb(abs(lb) < 1e-6) = 0;
0266 ub = Inf(nl,1);
0267 prob = update_Albub(prob, A, lb, ub);
0268
0269
0270 r = rxpp.';
0271 c = vars.xp.first - 1 + cxpp.';
0272 v = vxpp.';
0273 if ~nox
0274 r = [r, rxqp.'];
0275 c = [c, vars.xq.first - 1 + cxqp.'];
0276 v = [v, vxqp.'];
0277 end
0278 if branch_slack
0279 r = [r, 1:nl];
0280 c = [c, vars.sp.first:vars.sp.last];
0281 v = [v, -ones(1,nl)];
0282 end
0283 A = sparse(r, c, v, nl, total_vars);
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295 lb = -Inf(nl,1);
0296 ub = abs(Smax.*cos(alpha)) - Pl0;
0297 ub(abs(ub) < 1e-6) = 0;
0298 prob = update_Albub(prob, A, lb, ub);
0299
0300
0301 r = rxqq.';
0302 c = vars.xq.first - 1 + cxqq.';
0303 v = vxqq.';
0304 if ~nox
0305 r = [r, rxpq.'];
0306 c = [c, vars.xp.first - 1 + cxpq.'];
0307 v = [v, vxpq.'];
0308 end
0309 if branch_slack
0310 r = [r, 1:nl];
0311 c = [c, vars.sq.first:vars.sq.last];
0312 v = [v, ones(1,nl)];
0313 end
0314 A = sparse(r, c, v, nl, total_vars);
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326 lb = -abs(Smax.*sin(alpha)) - Ql0;
0327 lb(abs(lb) < 1e-6) = 0;
0328 ub = Inf(nl,1);
0329 prob = update_Albub(prob, A, lb, ub);
0330
0331
0332 r = rxqq.';
0333 c = vars.xq.first - 1 + cxqq.';
0334 v = vxqq.';
0335 if ~nox
0336 r = [r, rxpq.'];
0337 c = [c, vars.xp.first - 1 + cxpq.'];
0338 v = [v, vxpq.'];
0339 end
0340 if branch_slack
0341 r = [r, 1:nl];
0342 c = [c, vars.sq.first:vars.sq.last];
0343 v = [v, -ones(1,nl)];
0344 end
0345 A = sparse(r, c, v, nl, total_vars);
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357 lb = -Inf(nl,1);
0358 ub = abs(Smax.*sin(alpha)) - Ql0;
0359 ub(abs(ub) < 1e-6) = 0;
0360 prob = update_Albub(prob, A, lb, ub);
0361
0362
0363
0364
0365
0366 A = sparse(1,vars.xp.first:vars.xp.last, 1, 1, total_vars);
0367 lb = 0;
0368 ub = 0;
0369 prob = update_Albub(prob, A, lb, ub);
0370
0371 A = sparse(1,vars.xq.first:vars.xq.last, 1, 1, total_vars);
0372 lb = 0;
0373 ub = 0;
0374 prob = update_Albub(prob, A, lb, ub);
0375
0376
0377
0378
0379 A = sparse([1:2*nb,1:2*nb],...
0380 [vars.absx.first:vars.absx.last, vars.xp.first:vars.xp.last, vars.xq.first:vars.xq.last],...
0381 1, 2*nb, total_vars);
0382 lb = zeros(2*nb,1);
0383 ub = Inf(2*nb,1);
0384 prob = update_Albub(prob, A, lb, ub);
0385
0386 A = sparse([1:2*nb,1:2*nb],...
0387 [vars.absx.first:vars.absx.last, vars.xp.first:vars.xp.last, vars.xq.first:vars.xq.last], ...
0388 [ones(1,2*nb), -ones(1,2*nb)], 2*nb, total_vars);
0389 lb = zeros(2*nb,1);
0390 ub = Inf(2*nb,1);
0391 prob = update_Albub(prob, A, lb, ub);
0392
0393
0394 prob.xmin = -Inf(total_vars,1);
0395 prob.xmin(vars.xp.first:vars.xp.last) = -dPmax;
0396 prob.xmin(vars.xq.first:vars.xq.last) = -dQmax;
0397 prob.xmin(vars.absx.first:vars.absx.last) = 0;
0398 if branch_slack
0399 prob.xmin(vars.sp.first:vars.sp.last) = 0;
0400 prob.xmin(vars.sq.first:vars.sq.last) = 0;
0401 end
0402 if usedv
0403 prob.xmin(vars.sv.first:vars.sv.last) = 0;
0404 end
0405
0406 prob.xmax = Inf(total_vars,1);
0407 prob.xmax(vars.xp.first:vars.xp.last) = dPmax;
0408 prob.xmax(vars.xq.first:vars.xq.last) = dQmax;
0409
0410
0411
0412 r = vars.absx.first:vars.absx.last;
0413 v = ones(1, 2*nb);
0414 if branch_slack
0415 r = [r, vars.sp.first:vars.sp.last, vars.sq.first:vars.sq.last];
0416 v = [v, wspq*ones(1,2*nl)];
0417 end
0418 if usedv
0419 r = [r, vars.sv.first:vars.sv.last];
0420 v = [v, wsv*ones(1, vars.sv.last-vars.sv.first+1)];
0421 end
0422 prob.c = sparse(r, 1, v, total_vars, 1);
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441 prob.x0 = zeros(total_vars,1);
0442 if branch_slack
0443 prob.x0(vars.sp.first:vars.sp.last) = max(0, abs(Pl0) - abs(Smax.*cos(alpha)));
0444 prob.x0(vars.sq.first:vars.sq.last) = max(0, abs(Ql0) - abs(Smax.*sin(alpha)));
0445 end
0446 if usedv
0447 prob.x0(vars.sv.first) = max(0, max(max(Vmin - V0, V0 - Vmax)));
0448
0449 end
0450 if opt.vm.nodeperm.verbose > 1
0451 fprintf('\t problem setup time %0.3f sec\n', toc(tsetup));
0452 end
0453
0454 tsolve = tic;
0455 if opt.vm.nodeperm.verbose > 2
0456 prob.opt.verbose = opt.vm.nodeperm.verbose;
0457 end
0458 prob.opt.grb_opt.BarConvTol = 1e-4;
0459 [x, f, eflag, output, lambda] = qps_matpower(prob);
0460 if opt.vm.nodeperm.verbose > 1
0461 fprintf('\t solve time %0.3f sec\n', toc(tsolve));
0462 end
0463 if eflag
0464 dPb = x(vars.xp.first:vars.xp.last);
0465 dQb = x(vars.xq.first:vars.xq.last);
0466 Pbus = struct('old', real(Sbus), 'new', real(Sbus) + dPb);
0467 Qbus = struct('old', imag(Sbus), 'new', imag(Sbus) + dQb);
0468 else
0469 if ~branch_slack
0470 if opt.verbose > 0
0471 warning('sgvm_calc_injection_delta: optimization did not converge. Turning on branch slacks and retrying.')
0472 end
0473 opt.vm.nodeperm.branch_slack = true;
0474 [Pbus, Qbus] = sgvm_calc_injection_delta(mpc, opt);
0475 else
0476 error('sgvm_calc_injection_delta: optimization failed to converge.')
0477 end
0478 end
0479
0480 function prob = update_Albub(prob, A, lb, ub)
0481
0482
0483 prob.A = [prob.A; A];
0484 prob.l = [prob.l; lb];
0485 prob.u = [prob.u; ub];