0001 function mpc = toggle_dcline(mpc, on_off)
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 if strcmp(upper(on_off), 'ON')
0040
0041 c = idx_dcline;
0042
0043
0044 if ~isfield(mpc, 'dcline') || size(mpc.dcline, 2) < c.LOSS1
0045 error('toggle_dcline: case must contain a ''dcline'' field, an ndc x %d matrix.', c.LOSS1);
0046 end
0047 if isfield(mpc, 'dclinecost') && size(mpc.dcline, 1) ~= size(mpc.dclinecost, 1)
0048 error('toggle_dcline: number of rows in ''dcline'' field (%d) and ''dclinecost'' field (%d) do not match.', ...
0049 size(mpc.dcline, 1), size(mpc.dclinecost, 1));
0050 end
0051 l0 = mpc.dcline(:, c.LOSS0);
0052 l1 = mpc.dcline(:, c.LOSS1);
0053 k = find( (l0 + l1 .* mpc.dcline(:, c.PMIN) < 0) | ...
0054 (l0 + l1 .* mpc.dcline(:, c.PMAX) < 0) );
0055 if ~isempty(k)
0056 warning('toggle_dcline: loss can be negative for DC line from bus %d to %d\n', ...
0057 [mpc.dcline(k, c.F_BUS:c.T_BUS)]');
0058 end
0059
0060
0061
0062
0063 mpc = add_userfcn(mpc, 'ext2int', @userfcn_dcline_ext2int);
0064 mpc = add_userfcn(mpc, 'formulation', @userfcn_dcline_formulation);
0065 mpc = add_userfcn(mpc, 'int2ext', @userfcn_dcline_int2ext);
0066 mpc = add_userfcn(mpc, 'printpf', @userfcn_dcline_printpf);
0067 mpc = add_userfcn(mpc, 'savecase', @userfcn_dcline_savecase);
0068 mpc.userfcn.status.dcline = 1;
0069 elseif strcmp(upper(on_off), 'OFF')
0070 mpc = remove_userfcn(mpc, 'savecase', @userfcn_dcline_savecase);
0071 mpc = remove_userfcn(mpc, 'printpf', @userfcn_dcline_printpf);
0072 mpc = remove_userfcn(mpc, 'int2ext', @userfcn_dcline_int2ext);
0073 mpc = remove_userfcn(mpc, 'formulation', @userfcn_dcline_formulation);
0074 mpc = remove_userfcn(mpc, 'ext2int', @userfcn_dcline_ext2int);
0075 mpc.userfcn.status.dcline = 0;
0076 elseif strcmp(upper(on_off), 'STATUS')
0077 if isfield(mpc, 'userfcn') && isfield(mpc.userfcn, 'status') && ...
0078 isfield(mpc.userfcn.status, 'dcline')
0079 mpc = mpc.userfcn.status.dcline;
0080 else
0081 mpc = 0;
0082 end
0083 else
0084 error('toggle_dcline: 2nd argument must be ''on'', ''off'' or ''status''');
0085 end
0086
0087
0088
0089 function mpc = userfcn_dcline_ext2int(mpc, mpopt, args)
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0103 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0104 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0105 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0106 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0107 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0108 c = idx_dcline;
0109
0110
0111 if isfield(mpc, 'dclinecost')
0112 havecost = 1;
0113 else
0114 havecost = 0;
0115 end
0116
0117
0118 mpc.order.ext.dcline = mpc.dcline;
0119 if havecost
0120 mpc.order.ext.dclinecost = mpc.dclinecost;
0121 end
0122
0123
0124 mpc.order.dcline.status.on = find(mpc.dcline(:, c.BR_STATUS) > 0);
0125 mpc.order.dcline.status.off = find(mpc.dcline(:, c.BR_STATUS) <= 0);
0126
0127
0128 dc = mpc.dcline(mpc.order.dcline.status.on, :);
0129 if havecost
0130 dcc = mpc.dclinecost(mpc.order.dcline.status.on, :);
0131 mpc.dclinecost = dcc;
0132 end
0133 ndc = size(dc, 1);
0134 o = mpc.order;
0135
0136
0137 dc(:, c.F_BUS) = o.bus.e2i(dc(:, c.F_BUS));
0138 dc(:, c.T_BUS) = o.bus.e2i(dc(:, c.T_BUS));
0139 mpc.dcline = dc;
0140
0141
0142
0143
0144 dc(:, c.PT) = dc(:, c.PF) - (dc(:, c.LOSS0) + dc(:, c.LOSS1) .* dc(:, c.PF));
0145
0146
0147 fg = zeros(ndc, size(mpc.gen, 2));
0148 fg(:, MBASE) = 100;
0149 fg(:, GEN_STATUS) = dc(:, c.BR_STATUS);
0150 fg(:, PMIN) = -Inf;
0151 fg(:, PMAX) = Inf;
0152 tg = fg;
0153 fg(:, GEN_BUS) = dc(:, c.F_BUS);
0154 tg(:, GEN_BUS) = dc(:, c.T_BUS);
0155 fg(:, PG) = -dc(:, c.PF);
0156 tg(:, PG) = dc(:, c.PT);
0157 fg(:, QG) = dc(:, c.QF);
0158 tg(:, QG) = dc(:, c.QT);
0159 fg(:, VG) = dc(:, c.VF);
0160 tg(:, VG) = dc(:, c.VT);
0161 k = find(dc(:, c.PMIN) >= 0);
0162 if ~isempty(k)
0163 fg(k, PMAX) = -dc(k, c.PMIN);
0164 end
0165 k = find(dc(:, c.PMAX) >= 0);
0166 if ~isempty(k)
0167 fg(k, PMIN) = -dc(k, c.PMAX);
0168 end
0169 k = find(dc(:, c.PMIN) < 0);
0170 if ~isempty(k)
0171 tg(k, PMIN) = dc(k, c.PMIN);
0172 end
0173 k = find(dc(:, c.PMAX) < 0);
0174 if ~isempty(k)
0175 tg(k, PMAX) = dc(k, c.PMAX);
0176 end
0177 fg(:, QMIN) = dc(:, c.QMINF);
0178 fg(:, QMAX) = dc(:, c.QMAXF);
0179 tg(:, QMIN) = dc(:, c.QMINT);
0180 tg(:, QMAX) = dc(:, c.QMAXT);
0181
0182
0183
0184 fg(isload(fg), PMAX) = -1e-6;
0185 tg(isload(tg), PMAX) = -1e-6;
0186
0187
0188 refbus = find(mpc.bus(:, BUS_TYPE) == REF);
0189 mpc.bus(dc(:, c.F_BUS), BUS_TYPE) = PV;
0190 mpc.bus(dc(:, c.T_BUS), BUS_TYPE) = PV;
0191 mpc.bus(refbus, BUS_TYPE) = REF;
0192
0193
0194 nb = size(mpc.bus, 1);
0195 ng = size(mpc.gen, 1);
0196 if isfield(mpc, 'A') && ~isempty(mpc.A)
0197 [mA, nA] = size(mpc.A);
0198 if nA >= 2*nb + 2*ng
0199 mpc.A = [ mpc.A(:, 1:2*nb+ng) sparse(mA, 2*ndc) ...
0200 mpc.A(:, 2*nb+ng+(1:ng)) sparse(mA, 2*ndc) ...
0201 mpc.A(:, (2*nb+2*ng+1:nA)) ];
0202 else
0203 mpc.A = [ mpc.A(:, 1:nb+ng) sparse(mA, 2*ndc) ...
0204 mpc.A(:, (nb+ng+1:nA)) ];
0205 end
0206 end
0207 if isfield(mpc, 'N') && ~isempty(mpc.N)
0208 [mN, nN] = size(mpc.N);
0209 if nN >= 2*nb + 2*ng
0210 mpc.N = [ mpc.N(:, 1:2*nb+ng) sparse(mN, 2*ndc) ...
0211 mpc.N(:, 2*nb+ng+(1:ng)) sparse(mN, 2*ndc) ...
0212 mpc.N(:, (2*nb+2*ng+1:nN)) ];
0213 else
0214 mpc.N = [ mpc.N(:, 1:nb+ng) sparse(mN, 2*ndc) ...
0215 mpc.N(:, (nb+ng+1:nN)) ];
0216 end
0217 end
0218
0219
0220 mpc.gen = [mpc.gen; fg; tg];
0221
0222
0223 if isfield(mpc, 'gencost') && ~isempty(mpc.gencost)
0224 ngcc = size(mpc.gencost, 2);
0225 [pc, qc] = pqcost(mpc.gencost, ng);
0226 if havecost
0227 ndccc = size(dcc, 2);
0228 ccc = max([ngcc; ndccc]);
0229 if ccc > ngcc
0230 pc = [pc zeros(ng, ccc-ngcc)];
0231 if ~isempty(qc)
0232 qc = [qc zeros(ng, ccc-ngcc)];
0233 end
0234 ngcc = ccc;
0235 end
0236
0237
0238
0239 for k = 1:ndc
0240 if dcc(k, MODEL) == POLYNOMIAL
0241 nc = dcc(k, NCOST);
0242 temp = dcc(k, NCOST+(1:nc));
0243
0244
0245
0246 temp((nc-1):-2:1) = -temp((nc-1):-2:1);
0247 else
0248 nc = dcc(k, NCOST);
0249 temp = dcc(k, NCOST+(1:2*nc));
0250
0251 xx = -temp(1:2:2*nc);
0252 yy = temp(2:2:2*nc);
0253 temp(1:2:2*nc) = xx(end:-1:1);
0254 temp(2:2:2*nc) = yy(end:-1:1);
0255 end
0256 padding = zeros(1, ngcc-NCOST-length(temp));
0257 gck = [dcc(k, 1:NCOST) temp padding];
0258
0259
0260 pc = [pc; gck];
0261 end
0262
0263 dcgc = ones(ndc, 1) * [2 0 0 2 zeros(1, ngcc-4)];
0264 else
0265
0266 dcgc = ones(ndc, 1) * [2 0 0 2 zeros(1, ngcc-4)];
0267 pc = [pc; dcgc];
0268 end
0269
0270 pc = [pc; dcgc];
0271
0272
0273 if ~isempty(qc)
0274 qc = [qc; dcgc; dcgc];
0275 mpc.gencost = [pc; qc];
0276 else
0277 mpc.gencost = pc;
0278 end
0279 end
0280
0281
0282
0283 function om = userfcn_dcline_formulation(om, mpopt, args)
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309 c = idx_dcline;
0310
0311
0312 mpc = om.get_mpc();
0313 dc = mpc.dcline;
0314 ndc = size(dc, 1);
0315 ng = size(mpc.gen, 1) - 2*ndc;
0316
0317
0318 nL0 = -dc(:, c.LOSS0) / mpc.baseMVA;
0319 L1 = dc(:, c.LOSS1);
0320 Adc = [sparse(ndc, ng) spdiags(1-L1, 0, ndc, ndc) speye(ndc, ndc)];
0321
0322
0323 om.add_lin_constraint('dcline', Adc, nL0, nL0, {'Pg'});
0324
0325
0326
0327 function results = userfcn_dcline_int2ext(results, mpopt, args)
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0344 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0345 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0346 c = idx_dcline;
0347
0348
0349 o = results.order;
0350 k = find(o.ext.dcline(:, c.BR_STATUS));
0351 ndc = length(k);
0352 ng = size(results.gen, 1) - 2*ndc;
0353 nb = size(results.bus, 1);
0354
0355
0356 fg = results.gen(ng +(1:ndc), :);
0357 tg = results.gen(ng+ndc+(1:ndc), :);
0358
0359
0360 results.gen = results.gen(1:ng, :);
0361 if isfield(results, 'gencost') && ~isempty(results.gencost)
0362 results.gencost = results.gencost(1:ng, :);
0363 end
0364
0365
0366 if isfield(results, 'A') && ~isempty(results.A)
0367 [mA, nA] = size(results.A);
0368 if nA >= 2*nb + 2*ng + 4*ndc
0369 results.A = results.A(:, [1:2*nb+ng 2*nb+ng+2*ndc+(1:ng) 2*nb+2*ng+4*ndc+1:nA]);
0370 else
0371 results.A = results.A(:, [1:nb+ng nb+ng+2*ndc+1:nA]);
0372 end
0373 end
0374 if isfield(results, 'N') && ~isempty(results.N)
0375 [mN, nN] = size(results.N);
0376 if nN >= 2*nb + 2*ng + 4*ndc
0377 results.N = results.N(:, [1:2*nb+ng 2*nb+ng+2*ndc+(1:ng) 2*nb+2*ng+4*ndc+1:nN]);
0378 else
0379 results.N = results.N(:, [1:nb+ng nb+ng+2*ndc+1:nN]);
0380 end
0381 end
0382
0383
0384 results.dcline(:, c.PF) = -fg(:, PG);
0385 results.dcline(:, c.PT) = tg(:, PG);
0386 results.dcline(:, c.QF) = fg(:, QG);
0387 results.dcline(:, c.QT) = tg(:, QG);
0388 results.dcline(:, c.VF) = fg(:, VG);
0389 results.dcline(:, c.VT) = tg(:, VG);
0390 if size(fg, 2) >= MU_QMIN
0391 results.dcline(:, c.MU_PMIN ) = fg(:, MU_PMAX) + tg(:, MU_PMIN);
0392 results.dcline(:, c.MU_PMAX ) = fg(:, MU_PMIN) + tg(:, MU_PMAX);
0393 results.dcline(:, c.MU_QMINF) = fg(:, MU_QMIN);
0394 results.dcline(:, c.MU_QMAXF) = fg(:, MU_QMAX);
0395 results.dcline(:, c.MU_QMINT) = tg(:, MU_QMIN);
0396 results.dcline(:, c.MU_QMAXT) = tg(:, MU_QMAX);
0397 end
0398
0399
0400 results.order.int.dcline = results.dcline;
0401
0402 o.ext.dcline(k, c.PF:c.VT) = results.dcline(:, c.PF:c.VT);
0403 if size(results.dcline, 2) == c.MU_QMAXT
0404 o.ext.dcline(k, c.MU_PMIN:c.MU_QMAXT) = results.dcline(:, c.MU_PMIN:c.MU_QMAXT);
0405 end
0406 results.dcline = o.ext.dcline;
0407
0408
0409
0410 function results = userfcn_dcline_printpf(results, fd, mpopt, args)
0411
0412
0413
0414
0415
0416
0417
0418
0419 c = idx_dcline;
0420
0421
0422 SUPPRESS = mpopt.out.suppress_detail;
0423 if SUPPRESS == -1
0424 if size(results.bus, 1) > 500
0425 SUPPRESS = 1;
0426 else
0427 SUPPRESS = 0;
0428 end
0429 end
0430 OUT_ALL = mpopt.out.all;
0431 OUT_BRANCH = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.branch);
0432 if OUT_ALL == -1
0433 OUT_ALL_LIM = ~SUPPRESS * mpopt.out.lim.all;
0434 elseif OUT_ALL == 1
0435 OUT_ALL_LIM = 2;
0436 else
0437 OUT_ALL_LIM = 0;
0438 end
0439 if OUT_ALL_LIM == -1
0440 OUT_LINE_LIM = ~SUPPRESS * mpopt.out.lim.line;
0441 else
0442 OUT_LINE_LIM = OUT_ALL_LIM;
0443 end
0444 ctol = mpopt.opf.violation;
0445 ptol = 1e-4;
0446
0447
0448 dc = results.dcline;
0449 ndc = size(dc, 1);
0450 kk = find(dc(:, c.BR_STATUS) ~= 0);
0451 if OUT_BRANCH
0452 fprintf(fd, '\n================================================================================');
0453 fprintf(fd, '\n| DC Line Data |');
0454 fprintf(fd, '\n================================================================================');
0455 fprintf(fd, '\n Line From To Power Flow Loss Reactive Inj (MVAr)');
0456 fprintf(fd, '\n # Bus Bus From (MW) To (MW) (MW) From To ');
0457 fprintf(fd, '\n------ ------ ------ --------- --------- --------- --------- ---------');
0458 loss = 0;
0459 for k = 1:ndc
0460 if dc(k, c.BR_STATUS)
0461 fprintf(fd, '\n%5d%8d%8d%11.2f%11.2f%11.2f%11.2f%11.2f', ...
0462 k, dc(k, c.F_BUS:c.T_BUS), dc(k, c.PF:c.PT), ...
0463 dc(k, c.PF) - dc(k, c.PT), dc(k, c.QF:c.QT) );
0464 loss = loss + dc(k, c.PF) - dc(k, c.PT);
0465 else
0466 fprintf(fd, '\n%5d%8d%8d%11s%11s%11s%11s%11s', ...
0467 k, dc(k, c.F_BUS:c.T_BUS), '- ', '- ', '- ', '- ', '- ');
0468 end
0469 end
0470 fprintf(fd, '\n ---------');
0471 fprintf(fd, '\n Total:%11.2f\n', loss);
0472 end
0473
0474 if OUT_LINE_LIM == 2 || (OUT_LINE_LIM == 1 && ...
0475 (any(dc(kk, c.PF) > dc(kk, c.PMAX) - ctol) || ...
0476 any(dc(kk, c.MU_PMIN) > ptol) || ...
0477 any(dc(kk, c.MU_PMAX) > ptol)))
0478 fprintf(fd, '\n================================================================================');
0479 fprintf(fd, '\n| DC Line Constraints |');
0480 fprintf(fd, '\n================================================================================');
0481 fprintf(fd, '\n Line From To Minimum Actual Flow Maximum');
0482 fprintf(fd, '\n # Bus Bus Pmin mu Pmin (MW) Pmax Pmax mu ');
0483 fprintf(fd, '\n------ ------ ------ --------- --------- --------- --------- ---------');
0484 for k = 1:ndc
0485 if OUT_LINE_LIM == 2 || (OUT_LINE_LIM == 1 && ...
0486 (dc(k, c.PF) > dc(k, c.PMAX) - ctol || ...
0487 dc(k, c.MU_PMIN) > ptol || ...
0488 dc(k, c.MU_PMAX) > ptol))
0489 if dc(k, c.BR_STATUS)
0490 fprintf(fd, '\n%5d%8d%8d', k, dc(k, c.F_BUS:c.T_BUS) );
0491 if dc(k, c.MU_PMIN) > ptol
0492 fprintf(fd, '%11.3f', dc(k, c.MU_PMIN) );
0493 else
0494 fprintf(fd, '%11s', '- ' );
0495 end
0496 fprintf(fd, '%11.2f%11.2f%11.2f', ...
0497 dc(k, c.PMIN), dc(k, c.PF), dc(k, c.PMAX) );
0498 if dc(k, c.MU_PMAX) > ptol
0499 fprintf(fd, '%11.3f', dc(k, c.MU_PMAX) );
0500 else
0501 fprintf(fd, '%11s', '- ' );
0502 end
0503 else
0504 fprintf(fd, '\n%5d%8d%8d%11s%11s%11s%11s%11s', ...
0505 k, dc(k, c.F_BUS:c.T_BUS), '- ', '- ', '- ', '- ', '- ');
0506 end
0507 end
0508 end
0509 fprintf(fd, '\n');
0510 end
0511
0512
0513
0514 function mpc = userfcn_dcline_savecase(mpc, fd, prefix, args)
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524 c = idx_dcline;
0525
0526
0527 ncols = size(mpc.dcline, 2);
0528 fprintf(fd, '\n%%%%----- DC Line Data -----%%%%\n');
0529 if ncols < c.MU_QMAXT
0530 fprintf(fd, '%%\tfbus\ttbus\tstatus\tPf\tPt\tQf\tQt\tVf\tVt\tPmin\tPmax\tQminF\tQmaxF\tQminT\tQmaxT\tloss0\tloss1\n');
0531 else
0532 fprintf(fd, '%%\tfbus\ttbus\tstatus\tPf\tPt\tQf\tQt\tVf\tVt\tPmin\tPmax\tQminF\tQmaxF\tQminT\tQmaxT\tloss0\tloss1\tmuPmin\tmuPmax\tmuQminF\tmuQmaxF\tmuQminT\tmuQmaxT\n');
0533 end
0534 template = '\t%d\t%d\t%d\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g\t%.9g';
0535 if ncols == c.MU_QMAXT
0536 template = [template, '\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f'];
0537 end
0538 template = [template, ';\n'];
0539 fprintf(fd, '%sdcline = [\n', prefix);
0540 fprintf(fd, template, mpc.dcline.');
0541 fprintf(fd, '];\n');
0542
0543
0544