0001 function mpc = toggle_softlims(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
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
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107 if strcmp(upper(on_off), 'ON')
0108
0109
0110
0111
0112
0113
0114 mpc = add_userfcn(mpc, 'ext2int', @userfcn_softlims_ext2int);
0115 mpc = add_userfcn(mpc, 'formulation', @userfcn_softlims_formulation);
0116 mpc = add_userfcn(mpc, 'int2ext', @userfcn_softlims_int2ext);
0117 mpc = add_userfcn(mpc, 'printpf', @userfcn_softlims_printpf);
0118 mpc = add_userfcn(mpc, 'savecase', @userfcn_softlims_savecase);
0119 mpc.userfcn.status.softlims = 1;
0120 elseif strcmp(upper(on_off), 'OFF')
0121 mpc = remove_userfcn(mpc, 'savecase', @userfcn_softlims_savecase);
0122 mpc = remove_userfcn(mpc, 'printpf', @userfcn_softlims_printpf);
0123 mpc = remove_userfcn(mpc, 'int2ext', @userfcn_softlims_int2ext);
0124 mpc = remove_userfcn(mpc, 'formulation', @userfcn_softlims_formulation);
0125 mpc = remove_userfcn(mpc, 'ext2int', @userfcn_softlims_ext2int);
0126 mpc.userfcn.status.softlims = 0;
0127 elseif strcmp(upper(on_off), 'STATUS')
0128 if isfield(mpc, 'userfcn') && isfield(mpc.userfcn, 'status') && ...
0129 isfield(mpc.userfcn.status, 'softlims')
0130 mpc = mpc.userfcn.status.softlims;
0131 else
0132 mpc = 0;
0133 end
0134 else
0135 error('toggle_softlims: 2nd argument must be ''on'', ''off'' or ''status''');
0136 end
0137
0138
0139
0140 function mpc = userfcn_softlims_ext2int(mpc, mpopt, args)
0141
0142
0143
0144
0145
0146
0147
0148
0149 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0150 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0151 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0152 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0153 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0154 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0155 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0156 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0157 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0158
0159
0160 if ~isempty(mpopt) && isfield(mpopt, 'opf') && mpopt.opf.v_cartesian
0161 error('userfcn_softlims_ext2int: TOGGLE_SOFTLIMS is not implemented for OPF with cartesian voltages. Please set the MPOPT.opf.v_cartesian option to 0 for use with TOGGLE_SOFTLIMS.');
0162 end
0163
0164
0165 lims = softlims_lim2mat(mpopt);
0166 mat2lims = softlims_mat2lims(mpopt);
0167
0168
0169 mpc.softlims = softlims_defaults(mpc, mpopt);
0170
0171
0172 slims = softlims_init(mpc, mpopt);
0173 o = mpc.order;
0174
0175
0176 mpc.order.ext.softlims = slims;
0177
0178
0179 for m = fieldnames(mat2lims).'
0180 mat_name = m{:};
0181 n0 = size(o.ext.(mat_name), 1);
0182 n = size(mpc.(mat_name), 1);
0183 e2i = zeros(n0, 1);
0184 if strcmp(mat_name, 'gen')
0185
0186 e2i(o.gen.status.on(o.gen.i2e)) = (1:n)';
0187 else
0188 e2i(o.(mat_name).status.on) = (1:n)';
0189 end
0190 for lm = mat2lims.(mat_name)
0191 lim = lm{:};
0192 s = slims.(lim);
0193 if ~strcmp( s.hl_mod, 'none')
0194 s.idx = e2i(s.idx);
0195 k = find(s.idx == 0);
0196 s.idx(k) = [];
0197 s.cost(k, :) = [];
0198 s.sav(k) = [];
0199 s.ub(k) = [];
0200 if isfield(s, 'hl_val') && ~isscalar(s.hl_val)
0201 s.hl_val(k) = [];
0202 end
0203 slims.(lim) = s;
0204 end
0205 end
0206 end
0207
0208
0209
0210
0211 for lm = fieldnames(lims).'
0212 lim = lm{:};
0213 s = slims.(lim);
0214 if ~strcmp(s.hl_mod, 'none')
0215 mat_name = lims.(lim).mat_name;
0216 col = lims.(lim).col;
0217 mpc.(mat_name)(s.idx, col) = s.rval;
0218 end
0219 end
0220
0221 mpc.softlims = slims;
0222 mpc.order.int.softlims = slims;
0223
0224
0225
0226 function om = userfcn_softlims_formulation(om, mpopt, args)
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0237 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0238 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0239
0240
0241 mpc = om.get_mpc();
0242 [baseMVA, bus, branch] = deal(mpc.baseMVA, mpc.bus, mpc.branch);
0243 nb = size(bus, 1);
0244 ng = size(mpc.gen, 1);
0245
0246
0247 lims = softlims_lim2mat(mpopt);
0248
0249
0250 om.userdata.mpopt = mpopt;
0251
0252
0253 for lm = fieldnames(lims).'
0254 lim = lm{:};
0255 s = mpc.softlims.(lim);
0256 if strcmp(s.hl_mod, 'none')
0257 continue;
0258 end
0259 varname = ['s_', lower(lim)];
0260 cstname = ['cs_', lower(lim)];
0261 ns = length(s.idx);
0262
0263
0264 switch lim
0265 case {'VMIN', 'VMAX'}
0266 ub = s.ub;
0267 Cw = s.cost(:, 1);
0268 case {'RATE_A', 'PMIN', 'PMAX', 'QMIN', 'QMAX'}
0269 ub = s.ub / baseMVA;
0270 Cw = s.cost(:, 1) * baseMVA;
0271 case {'ANGMIN', 'ANGMAX'}
0272 ub = s.ub * pi/180;
0273 Cw = s.cost(:, 1) * 180/pi;
0274 otherwise
0275 error('userfcn_soflims_formulation: limit %s is unknown ', lim)
0276 end
0277 om.add_var(varname, ns, zeros(ns, 1), zeros(ns, 1), ub);
0278 om.add_quad_cost(cstname, [], Cw, 0, {varname});
0279
0280
0281 switch lim
0282 case 'ANGMIN'
0283
0284 ns = length(s.idx);
0285 Av = sparse([1:ns,1:ns].', ...
0286 [branch(s.idx, F_BUS); branch(s.idx, T_BUS)], ...
0287 [ones(ns,1);-ones(ns,1)], ns, nb);
0288 As = speye(ns);
0289 lb = s.sav * pi/180;
0290 ub = Inf(ns,1);
0291
0292 om.add_lin_constraint('soft_angmin', [Av As], lb, ub, ...
0293 {'Va', 's_angmin'});
0294 case 'ANGMAX'
0295
0296 ns = length(s.idx);
0297 Av = sparse([1:ns,1:ns].', ...
0298 [branch(s.idx, F_BUS); branch(s.idx, T_BUS)], ...
0299 [ones(ns,1);-ones(ns,1)], ns, nb);
0300 As = speye(ns);
0301 lb = -Inf(ns,1);
0302 ub = s.sav * pi/180;
0303
0304 om.add_lin_constraint('soft_angmax', [Av -As], lb, ub, ...
0305 {'Va', 's_angmax'});
0306 case 'PMIN'
0307
0308 ns = length(s.idx);
0309 Av = sparse(1:ns, s.idx, 1, ns, ng);
0310 As = speye(ns);
0311 lb = s.sav / baseMVA;
0312 ub = Inf(ns,1);
0313
0314 om.add_lin_constraint('soft_pmin', [Av As], lb, ub, ...
0315 {'Pg', 's_pmin'});
0316 case 'PMAX'
0317
0318 ns = length(s.idx);
0319 Av = sparse(1:ns, s.idx, 1, ns, ng);
0320 As = speye(ns);
0321 lb = -Inf(ns,1);
0322 ub = s.sav / baseMVA;
0323
0324 om.add_lin_constraint('soft_pmax', [Av -As], lb, ub, ...
0325 {'Pg', 's_pmax'});
0326 end
0327 end
0328
0329
0330 isDC = strcmp(mpopt.model, 'DC');
0331 lims = {'RATE_A'};
0332 if ~isDC
0333 lims = {'VMIN', 'VMAX', 'QMIN', 'QMAX', lims{:}};
0334 end
0335 for lm = lims
0336 lim = lm{:};
0337 s = mpc.softlims.(lim);
0338 if strcmp(s.hl_mod, 'none')
0339 continue;
0340 end
0341 ns = length(s.idx);
0342
0343 switch lim
0344 case 'RATE_A'
0345 if isDC
0346
0347 Bf = om.get_userdata('Bf');
0348 Pfinj = om.get_userdata('Pfinj');
0349
0350
0351
0352
0353 I = speye(ns, ns);
0354 Asf = [ Bf(s.idx, :) -I];
0355 Ast = [-Bf(s.idx, :) -I];
0356 lsf = -Inf(ns, 1);
0357 lst = lsf;
0358 usf = -Pfinj(s.idx) + s.sav/baseMVA;
0359 ust = Pfinj(s.idx) + s.sav/baseMVA;
0360 vs = {'Va', 's_rate_a'};
0361
0362 om.add_lin_constraint('softPf', Asf, lsf, usf, vs);
0363 om.add_lin_constraint('softPt', Ast, lst, ust, vs);
0364 else
0365 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0366
0367 fcn = @(x)softlims_fcn(x, mpc, Yf(s.idx, :), Yt(s.idx, :), ...
0368 s.idx, mpopt, s.sav/baseMVA);
0369 hess = @(x, lam)softlims_hess(x, lam, mpc, Yf(s.idx, :), ...
0370 Yt(s.idx, :), s.idx, mpopt);
0371 om.add_nln_constraint({'softSf', 'softSt'}, [ns;ns], 0, ...
0372 fcn, hess, {'Va', 'Vm', 's_rate_a'});
0373 end
0374 case 'VMIN'
0375
0376 Av = sparse(1:ns, s.idx, 1, ns, nb);
0377 As = speye(ns);
0378 lb = s.sav;
0379 ub = Inf(ns, 1);
0380
0381 om.add_lin_constraint('soft_vmin', [Av As], lb, ub, ...
0382 {'Vm', 's_vmin'});
0383 case 'VMAX'
0384
0385 Av = sparse(1:ns, s.idx, 1, ns, nb);
0386 As = speye(ns);
0387 lb = -Inf(ns, 1);
0388 ub = s.sav;
0389
0390 om.add_lin_constraint('soft_vmax', [Av -As], lb, ub, ...
0391 {'Vm', 's_vmax'});
0392 case 'QMIN'
0393
0394 Av = sparse(1:ns, s.idx, 1, ns, ng);
0395 As = speye(ns);
0396 lb = s.sav / baseMVA;
0397 ub = Inf(ns, 1);
0398
0399 om.add_lin_constraint('soft_qmin', [Av As], lb, ub, ...
0400 {'Qg', 's_qmin'});
0401 case 'QMAX'
0402
0403 Av = sparse(1:ns, s.idx, 1, ns, ng);
0404 As = speye(ns);
0405 lb = -Inf(ns, 1);
0406 ub = s.sav / baseMVA;
0407
0408 om.add_lin_constraint('soft_qmax', [Av -As], lb, ub, ...
0409 {'Qg', 's_qmax'});
0410 end
0411 end
0412
0413
0414
0415 function results = userfcn_softlims_int2ext(results, mpopt, args)
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0439 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0440 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0441 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0442 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0443 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0444 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0445 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0446
0447
0448 slims = results.softlims;
0449 isOPF = isfield(results, 'f') && ~isempty(results.f);
0450 if isOPF
0451 mpopt = results.om.get_userdata('mpopt');
0452 results.om.userdata = rmfield(results.om.userdata, 'mpopt');
0453 isDC = strcmp(mpopt.model, 'DC');
0454 end
0455
0456
0457 lims = softlims_lim2mat(mpopt);
0458
0459
0460 o = results.order;
0461 results.softlims = o.ext.softlims;
0462
0463
0464 for lm = fieldnames(lims).'
0465 lim = lm{:};
0466 s = slims.(lim);
0467 if strcmp(s.hl_mod, 'none')
0468 continue;
0469 end
0470 mat_name = lims.(lim).mat_name;
0471 col = lims.(lim).col;
0472 results.(mat_name)(s.idx, col) = s.sav;
0473 end
0474
0475
0476 for lm = fieldnames(lims).'
0477 lim = lm{:};
0478 if strcmp(slims.(lim).hl_mod, 'none')
0479 continue;
0480 end
0481 results.softlims.(lim) = rmfield(results.softlims.(lim), 'sav');
0482 results.softlims.(lim) = rmfield(results.softlims.(lim), 'rval');
0483 results.softlims.(lim) = rmfield(results.softlims.(lim), 'ub');
0484 end
0485
0486
0487
0488 tol = 1e-8;
0489 if isOPF
0490 for lm = fieldnames(lims).'
0491 lim = lm{:};
0492 s = slims.(lim);
0493 if strcmp(s.hl_mod, 'none')
0494 continue;
0495 end
0496
0497 mat_name = lims.(lim).mat_name;
0498 n0 = size(o.ext.(mat_name), 1);
0499 n = size(results.(mat_name), 1);
0500 results.softlims.(lim).overload = zeros(n0, 1);
0501 results.softlims.(lim).ovl_cost = zeros(n0, 1);
0502 varname = ['s_', lower(lim)];
0503 if isfield(results.var.val, varname)
0504 var = results.var.val.(varname);
0505 var(var < tol) = 0;
0506 end
0507 switch lim
0508 case {'VMIN', 'VMAX'}
0509 if isDC
0510 continue;
0511 end
0512 case {'RATE_A', 'PMIN', 'PMAX'}
0513 var = var * results.baseMVA;
0514 case {'QMIN', 'QMAX'}
0515 if isDC
0516 continue;
0517 end
0518 var = var * results.baseMVA;
0519 case {'ANGMIN', 'ANGMAX'}
0520 var = var * 180/pi;
0521 otherwise
0522 error('userfcn_soflims_formulation: limit %s is unknown ', lim)
0523 end
0524
0525
0526
0527 if strcmp(mat_name, 'gen')
0528 results.softlims.(lim).overload(o.gen.status.on(o.gen.i2e(s.idx))) = var;
0529 results.softlims.(lim).ovl_cost(o.gen.status.on(o.gen.i2e(s.idx))) = var .* s.cost(:, 1);
0530 else
0531 results.softlims.(lim).overload(o.(mat_name).status.on(s.idx)) = var;
0532 results.softlims.(lim).ovl_cost(o.(mat_name).status.on(s.idx)) = var .* s.cost(:, 1);
0533 end
0534
0535 cname = ['soft_' lower(lim)];
0536 switch lim
0537 case {'ANGMAX', 'PMAX', 'VMAX', 'QMAX'}
0538 mu = results.lin.mu.u.(cname);
0539 case {'ANGMIN', 'PMIN', 'VMIN', 'QMIN'}
0540 mu = results.lin.mu.l.(cname);
0541 case 'RATE_A'
0542 if isDC
0543 muF = results.lin.mu.u.softPf;
0544 muT = results.lin.mu.u.softPt;
0545 else
0546 muF = results.nli.mu.softSf;
0547 muT = results.nli.mu.softSt;
0548 end
0549 end
0550 switch lim
0551 case 'ANGMAX'
0552 results.branch(slims.ANGMAX.idx, MU_ANGMAX) = mu * pi/180;
0553 case 'ANGMIN'
0554 results.branch(slims.ANGMIN.idx, MU_ANGMIN) = mu * pi/180;
0555 case 'PMAX'
0556 results.gen(slims.PMAX.idx, MU_PMAX) = mu / results.baseMVA;
0557 case 'PMIN'
0558 results.gen(slims.PMIN.idx, MU_PMIN) = mu / results.baseMVA;
0559 case 'RATE_A'
0560 results.branch(slims.RATE_A.idx, MU_SF) = muF / results.baseMVA;
0561 results.branch(slims.RATE_A.idx, MU_ST) = muT / results.baseMVA;
0562 end
0563 if ~isDC
0564 switch lim
0565 case 'RATE_A'
0566 if upper(mpopt.opf.flow_lim(1)) ~= 'P'
0567 s = slims.RATE_A;
0568 var = results.softlims.RATE_A.overload(o.branch.status.on(s.idx));
0569
0570 cf = 2 * (s.sav + var) / results.baseMVA;
0571 results.branch(s.idx, MU_SF) = ...
0572 results.branch(s.idx, MU_SF) .* cf;
0573 results.branch(s.idx, MU_ST) = ...
0574 results.branch(s.idx, MU_ST) .* cf;
0575 end
0576 case 'VMAX'
0577 results.bus(slims.VMAX.idx, MU_VMAX) = mu;
0578 case 'VMIN'
0579 results.bus(slims.VMIN.idx, MU_VMIN) = mu;
0580 case 'QMAX'
0581 results.gen(slims.QMAX.idx, MU_QMAX) = mu / results.baseMVA;
0582 case 'QMIN'
0583 results.gen(slims.QMIN.idx, MU_QMIN) = mu / results.baseMVA;
0584 end
0585 end
0586 end
0587 end
0588
0589
0590
0591 function results = userfcn_softlims_printpf(results, fd, mpopt, args)
0592
0593
0594
0595
0596
0597
0598
0599
0600 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0601 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0602 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0603 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0604 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0605 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0606 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0607 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0608
0609
0610 ptol = 1e-6;
0611 isOPF = isfield(results, 'f') && ~isempty(results.f);
0612 SUPPRESS = mpopt.out.suppress_detail;
0613 if SUPPRESS == -1
0614 if size(results.bus, 1) > 500
0615 SUPPRESS = 1;
0616 else
0617 SUPPRESS = 0;
0618 end
0619 end
0620 OUT_ALL = mpopt.out.all;
0621 OUT_FORCE = mpopt.out.force;
0622 OUT_V_LIM = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.v);
0623 OUT_LINE_LIM = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.line);
0624 OUT_PG_LIM = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.pg);
0625 OUT_QG_LIM = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.qg);
0626
0627 if isOPF && (results.success || OUT_FORCE)
0628 [bus, gen, branch] = deal(results.bus, results.gen, results.branch);
0629 isAC = strcmp(mpopt.model, 'AC');
0630 if isAC
0631 s = results.softlims.VMAX;
0632 if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0633 fprintf(fd, '\n================================================================================');
0634 fprintf(fd, '\n| Soft Voltage Upper Bounds |');
0635 fprintf(fd, '\n================================================================================');
0636 k = find(s.overload(s.idx) | bus(s.idx, MU_VMAX) > ptol);
0637 if isempty(k)
0638 fprintf(fd,'\nNo violations.\n');
0639 else
0640 fprintf(fd, '\nBus Voltage Limit Overload mu');
0641 fprintf(fd, '\n # Mag(pu) (pu) (pu) ($/pu)');
0642 fprintf(fd, '\n----- ------- ------- ------- ---------');
0643 fprintf(fd, '\n%5d%8.3f%9.3f%9.3f%11.3f',...
0644 [ bus(s.idx(k), BUS_I), bus(s.idx(k),[VM, VMAX]),...
0645 s.overload(s.idx(k)), ...
0646 bus(s.idx(k), MU_VMAX)...
0647 ]');
0648 fprintf(fd, '\n --------');
0649 fprintf(fd, '\n Total:%10.2f', ...
0650 sum(s.overload(s.idx(k))));
0651 fprintf(fd, '\n');
0652 end
0653 end
0654 s = results.softlims.VMIN;
0655 if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0656 fprintf(fd, '\n================================================================================');
0657 fprintf(fd, '\n| Soft Voltage Lower Bounds |');
0658 fprintf(fd, '\n================================================================================');
0659 k = find(s.overload(s.idx) | bus(s.idx, MU_VMIN) > ptol);
0660 if isempty(k)
0661 fprintf(fd,'\nNo violations.\n');
0662 else
0663 fprintf(fd, '\n----------------------------------------');
0664 fprintf(fd, '\n Bus Voltage Limit Overload mu');
0665 fprintf(fd, '\n # Mag(pu) (pu) (pu) ($/pu)');
0666 fprintf(fd, '\n----- ------- ------- ------- ---------');
0667 fprintf(fd, '\n%5d%8.3f%9.3f%9.3f%11.3f',...
0668 [ bus(s.idx(k), BUS_I), bus(s.idx(k),[VM, VMIN]),...
0669 s.overload(s.idx(k)), ...
0670 bus(s.idx(k), MU_VMIN)...
0671 ]');
0672 fprintf(fd, '\n --------');
0673 fprintf(fd, '\n Total:%10.2f', ...
0674 sum(s.overload(s.idx(k))));
0675 fprintf(fd, '\n');
0676 end
0677 end
0678 end
0679 s = results.softlims.PMAX;
0680 if OUT_PG_LIM && ~strcmp(s.hl_mod,'none')
0681 k = find(s.overload(s.idx) | gen(s.idx, MU_PMAX) > ptol);
0682 fprintf(fd, '\n================================================================================');
0683 fprintf(fd, '\n| Soft Generator Active Power Upper Bounds |');
0684 fprintf(fd, '\n================================================================================');
0685 if isempty(k)
0686 fprintf(fd,'\nNo violations.\n');
0687 else
0688 fprintf(fd, '\nGen Bus Generation Limit Overload mu');
0689 fprintf(fd, '\n # # P (MW) (MW) (MW) ($/MW)');
0690 fprintf(fd, '\n----- ----- -------- ------- ------- ---------');
0691 fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0692 [ s.idx(k), gen(s.idx(k),[GEN_BUS, PG, PMAX]),...
0693 s.overload(s.idx(k)), ...
0694 gen(s.idx(k), MU_PMAX)...
0695 ]');
0696 fprintf(fd, '\n --------');
0697 fprintf(fd, '\n Total:%10.2f', ...
0698 sum(s.overload(s.idx(k))));
0699 fprintf(fd, '\n');
0700 end
0701 end
0702 s = results.softlims.PMIN;
0703 if OUT_PG_LIM && ~strcmp(s.hl_mod,'none')
0704 k = find(s.overload(s.idx) | gen(s.idx, MU_PMIN) > ptol);
0705 fprintf(fd, '\n================================================================================');
0706 fprintf(fd, '\n| Soft Generator Active Power Lower Bounds |');
0707 fprintf(fd, '\n================================================================================');
0708 if isempty(k)
0709 fprintf(fd,'\nNo violations.\n');
0710 else
0711 fprintf(fd, '\nGen Bus Generation Limit Overload mu');
0712 fprintf(fd, '\n # # P (MW) (MW) (MW) ($/MW)');
0713 fprintf(fd, '\n----- ----- -------- ------- ------- ---------');
0714 fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0715 [ s.idx(k), gen(s.idx(k),[GEN_BUS, PG, PMIN]),...
0716 s.overload(s.idx(k)), ...
0717 gen(s.idx(k), MU_PMIN)...
0718 ]');
0719 fprintf(fd, '\n --------');
0720 fprintf(fd, '\n Total:%10.2f', ...
0721 sum(s.overload(s.idx(k))));
0722 fprintf(fd, '\n');
0723 end
0724 end
0725 if isAC
0726 s = results.softlims.QMAX;
0727 if OUT_QG_LIM && ~strcmp(s.hl_mod,'none')
0728 k = find(s.overload(s.idx) | gen(s.idx, MU_QMAX) > ptol);
0729 fprintf(fd, '\n================================================================================');
0730 fprintf(fd, '\n| Soft Generator Reactive Power Upper Bounds |');
0731 fprintf(fd, '\n================================================================================');
0732 if isempty(k)
0733 fprintf(fd,'\nNo violations.\n');
0734 else
0735 fprintf(fd, '\nGen Bus Generation Limit Overload mu');
0736 fprintf(fd, '\n # # Q (MVAr) Q (MVAr) (MVAr) ($/MVAr)');
0737 fprintf(fd, '\n----- ----- -------- ------- ------- ---------');
0738 fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0739 [ s.idx(k), gen(s.idx(k),[GEN_BUS, QG, QMAX]),...
0740 s.overload(s.idx(k)), ...
0741 gen(s.idx(k), MU_QMAX)...
0742 ]');
0743 fprintf(fd, '\n --------');
0744 fprintf(fd, '\n Total:%10.2f', ...
0745 sum(s.overload(s.idx(k))));
0746 fprintf(fd, '\n');
0747 end
0748 end
0749 s = results.softlims.QMIN;
0750 if OUT_QG_LIM && ~strcmp(s.hl_mod,'none')
0751 k = find(s.overload(s.idx) | gen(s.idx, MU_QMIN) > ptol);
0752 fprintf(fd, '\n================================================================================');
0753 fprintf(fd, '\n| Soft Generator Reactive Power Lower Bounds |');
0754 fprintf(fd, '\n================================================================================');
0755 if isempty(k)
0756 fprintf(fd,'\nNo violations.\n');
0757 else
0758 fprintf(fd, '\nGen Bus Generation Limit Overload mu');
0759 fprintf(fd, '\n # # Q (MVAr) Q (MVAr) (MVAr) ($/MVAr)');
0760 fprintf(fd, '\n----- ----- -------- ------- ------- ---------');
0761 fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0762 [ s.idx(k), gen(s.idx(k),[GEN_BUS, QG, QMIN]),...
0763 s.overload(s.idx(k)), ...
0764 gen(s.idx(k), MU_QMIN)...
0765 ]');
0766 fprintf(fd, '\n --------');
0767 fprintf(fd, '\n Total:%10.2f', ...
0768 sum(s.overload(s.idx(k))));
0769 fprintf(fd, '\n');
0770 end
0771 end
0772 end
0773 s = results.softlims.RATE_A;
0774 if OUT_LINE_LIM && ~strcmp(s.hl_mod, 'none')
0775 fprintf(fd, '\n================================================================================');
0776 fprintf(fd, '\n| Soft Branch Flow Limits |');
0777 fprintf(fd, '\n================================================================================');
0778 k = find(s.overload(s.idx) | sum(branch(s.idx, MU_SF:MU_ST), 2) > ptol);
0779 if isempty(k)
0780 fprintf(fd,'\nNo violations.\n');
0781 else
0782
0783 lim_type = upper(mpopt.opf.flow_lim(1));
0784 if ~isAC || lim_type == 'P' || lim_type == '2'
0785 Ff = branch(s.idx(k), PF);
0786 Ft = branch(s.idx(k), PT);
0787
0788 str = '\n # Bus Bus (MW) (MW) (MW) ($/MW)';
0789 elseif lim_type == 'I'
0790
0791 i2e = bus(:, BUS_I);
0792 e2i = sparse(max(i2e), 1);
0793 e2i(i2e) = (1:size(bus, 1))';
0794
0795 V = bus(:, VM) .* exp(1j * pi/180 * bus(:, VA));
0796 Ff = abs( (branch(s.idx(k), PF) + 1j * branch(s.idx(k), QF)) ./ V(e2i(branch(s.idx(k), F_BUS))) );
0797 Ft = abs( (branch(s.idx(k), PT) + 1j * branch(s.idx(k), QT)) ./ V(e2i(branch(s.idx(k), T_BUS))) );
0798
0799 str = '\n # Bus Bus (kA*basekV) (kA*basekV) (kA*basekV) ($/kA/basekV)';
0800 else
0801 Ff = abs(branch(s.idx(k), PF) + 1j * branch(s.idx(k), QF));
0802 Ft = abs(branch(s.idx(k), PT) + 1j * branch(s.idx(k), QT));
0803
0804 str = '\n # Bus Bus (MVA) (MVA) (MVA) ($/MVA)';
0805 end
0806 flow = max(abs(Ff), abs(Ft));
0807
0808
0809 fprintf(fd, '\nBrnch From To Flow Limit Overload mu');
0810 fprintf(fd, str);
0811
0812 fprintf(fd, '\n----- ----- ----- ----------- ----------- ----------- -------------');
0813
0814 fprintf(fd, '\n%4d%7d%7d%13.2f%13.2f%13.2f%14.3f', ...
0815 [ s.idx(k), branch(s.idx(k), [F_BUS, T_BUS]), ...
0816 flow, branch(s.idx(k), RATE_A), ...
0817 s.overload(s.idx(k)), ...
0818 sum(branch(s.idx(k), MU_SF:MU_ST), 2) ...
0819 ]');
0820 fprintf(fd, '\n -----------');
0821 fprintf(fd, '\n Total: %10.2f', ...
0822 sum(s.overload(s.idx(k))));
0823 fprintf(fd, '\n');
0824 end
0825 end
0826 s = results.softlims.ANGMAX;
0827 if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0828 k = find(s.overload(s.idx) | branch(s.idx, MU_ANGMAX) > ptol);
0829 delta = calc_branch_angle(results);
0830 fprintf(fd, '\n================================================================================');
0831 fprintf(fd, '\n| Soft Maximum Angle Difference Limits |');
0832 fprintf(fd, '\n================================================================================');
0833 if isempty(k)
0834 fprintf(fd,'\nNo violations.\n');
0835 else
0836 fprintf(fd, '\nBrnch From To Angle Limit Overload mu');
0837 fprintf(fd, '\n # Bus Bus (deg) (deg) (deg) ($/MW)');
0838 fprintf(fd, '\n----- ----- ----- -------- -------- -------- ---------');
0839 fprintf(fd, '\n%4d%7d%7d%10.3f%10.3f%10.3f%11.3f', ...
0840 [ s.idx(k), branch(s.idx(k), [F_BUS, T_BUS]), ...
0841 delta(s.idx(k)), branch(s.idx(k), ANGMAX),...
0842 s.overload(s.idx(k)),...
0843 branch(s.idx(k), MU_ANGMAX)...
0844 ]');
0845 fprintf(fd, '\n --------');
0846 fprintf(fd, '\n Total:%10.2f', ...
0847 sum(s.overload(s.idx(k))));
0848 fprintf(fd, '\n');
0849 end
0850 end
0851 s = results.softlims.ANGMIN;
0852 if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0853 k = find(s.overload(s.idx) | branch(s.idx, MU_ANGMIN) > ptol);
0854 delta = calc_branch_angle(results);
0855 fprintf(fd, '\n================================================================================');
0856 fprintf(fd, '\n| Soft Minimum Angle Difference Limits |');
0857 fprintf(fd, '\n================================================================================');
0858 if isempty(k)
0859 fprintf(fd,'\nNo violations.\n');
0860 else
0861 fprintf(fd, '\nBrnch From To Angle Limit Overload mu');
0862 fprintf(fd, '\n # Bus Bus (deg) (deg) (deg) ($/MW)');
0863 fprintf(fd, '\n----- ----- ----- -------- -------- -------- ---------');
0864 fprintf(fd, '\n%4d%7d%7d%10.3f%10.3f%10.3f%11.3f', ...
0865 [ s.idx(k), branch(s.idx(k), [F_BUS, T_BUS]), ...
0866 delta(s.idx(k)), branch(s.idx(k), ANGMIN),...
0867 s.overload(s.idx(k)),...
0868 branch(s.idx(k), MU_ANGMIN)...
0869 ]');
0870 fprintf(fd, '\n --------');
0871 fprintf(fd, '\n Total:%10.2f', ...
0872 sum(s.overload(s.idx(k))));
0873 fprintf(fd, '\n');
0874 end
0875 end
0876 end
0877
0878
0879
0880 function mpc = userfcn_softlims_savecase(mpc, fd, prefix, args)
0881
0882
0883
0884
0885
0886
0887
0888
0889
0890 lims = softlims_lim2mat();
0891
0892
0893 slfieldnames = { 'hl_mod', 'hl_val', 'idx', 'busnum', 'cost', ...
0894 'overload', 'ovl_cost' };
0895 fields = struct( ...
0896 'hl_mod', struct('desc', 'type of hard limit modification', 'tok', '%s'),...
0897 'hl_val', struct('desc', 'value(s) used to set new hard limit', 'tok', '%g'),...
0898 'busnum', struct('desc', 'bus numbers for soft voltage limit', 'tok', '%d'),...
0899 'idx', struct('desc', '%s matrix row indices', 'tok', '%d'),...
0900 'cost', struct('desc', 'violation cost coefficient', 'tok', '%g'),...
0901 'overload', struct('desc', 'violation of original limit', 'tok', '%0.10g'),...
0902 'ovl_cost', struct('desc', 'overload cost', 'tok', '%.10g') ...
0903 );
0904
0905 if isfield(mpc, 'softlims')
0906 fprintf(fd, '\n%%%%----- OPF Soft Limit Data -----%%%%\n');
0907 limits = fieldnames(lims).';
0908 for lm = limits
0909 lim = lm{:};
0910 s = mpc.softlims.(lim);
0911 if ~strcmp(lim, limits{1})
0912 fprintf(fd, '\n');
0913 end
0914 fprintf(fd,'%%%% %s soft limit data\n', lim);
0915 for f = slfieldnames
0916 f = f{1};
0917 if isfield(s, f)
0918 if strcmp(f, 'idx')
0919 if ismember(lim, {'VMIN', 'VMAX'})
0920 desc = fields.busnum.desc;
0921 else
0922 desc = sprintf(fields.idx.desc, lims.(lim).mat_name);
0923 end
0924 else
0925 desc = fields.(f).desc;
0926 end
0927 if strcmp(f, 'hl_mod')
0928 str = sprintf('%ssoftlims.%s.hl_mod = ''%s'';', prefix, lim, s.hl_mod);
0929 pad = repmat(' ', 1, 40-length(str));
0930 fprintf(fd, '%s%s%%%% %s\n', str, pad, desc);
0931 else
0932 str = sprintf('%ssoftlims.%s.%s = [', prefix, lim, f);
0933 pad = repmat(' ', 1, 40-length(str));
0934 fprintf(fd, '%s%s%%%% %s\n', str, pad, desc);
0935 fprintf(fd, ['\t', fields.(f).tok, ';\n'], s.(f));
0936 fprintf(fd, '];\n');
0937 end
0938 end
0939 end
0940 end
0941 end
0942
0943
0944
0945 function lim2mat = softlims_lim2mat(mpopt)
0946
0947
0948 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0949 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0950 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0951 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0952 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0953 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0954 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0955 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0956
0957 if nargin < 1
0958 mpopt = [];
0959 end
0960 if ~isempty(mpopt) && isfield(mpopt, 'model') && strcmp(mpopt.model, 'DC')
0961 lim2mat = struct(...
0962 'PMAX', struct('mat_name', 'gen', 'col', PMAX, 'dir', 1 ), ...
0963 'PMIN', struct('mat_name', 'gen', 'col', PMIN, 'dir', -1 ), ...
0964 'RATE_A', struct('mat_name', 'branch', 'col', RATE_A, 'dir', 1 ), ...
0965 'ANGMAX', struct('mat_name', 'branch', 'col', ANGMAX, 'dir', 1 ), ...
0966 'ANGMIN', struct('mat_name', 'branch', 'col', ANGMIN, 'dir', -1 ) ...
0967 );
0968 else
0969 lim2mat = struct(...
0970 'VMAX', struct('mat_name', 'bus', 'col', VMAX, 'dir', 1 ), ...
0971 'VMIN', struct('mat_name', 'bus', 'col', VMIN, 'dir', -1 ), ...
0972 'PMAX', struct('mat_name', 'gen', 'col', PMAX, 'dir', 1 ), ...
0973 'PMIN', struct('mat_name', 'gen', 'col', PMIN, 'dir', -1 ), ...
0974 'QMAX', struct('mat_name', 'gen', 'col', QMAX, 'dir', 1 ), ...
0975 'QMIN', struct('mat_name', 'gen', 'col', QMIN, 'dir', -1 ), ...
0976 'RATE_A', struct('mat_name', 'branch', 'col', RATE_A, 'dir', 1 ), ...
0977 'ANGMAX', struct('mat_name', 'branch', 'col', ANGMAX, 'dir', 1 ), ...
0978 'ANGMIN', struct('mat_name', 'branch', 'col', ANGMIN, 'dir', -1 ) ...
0979 );
0980 end
0981
0982
0983
0984 function mat2lims = softlims_mat2lims(mpopt)
0985 if nargin < 1
0986 mpopt = [];
0987 end
0988 if ~isempty(mpopt) && isfield(mpopt, 'model') && strcmp(mpopt.model, 'DC')
0989 mat2lims = struct(...
0990 'branch', {{'ANGMAX', 'ANGMIN', 'RATE_A'}}, ...
0991 'gen', {{'PMAX', 'PMIN'}} ...
0992 );
0993 else
0994 mat2lims = struct(...
0995 'bus', {{'VMAX', 'VMIN'}}, ...
0996 'branch', {{'ANGMAX', 'ANGMIN', 'RATE_A'}}, ...
0997 'gen', {{'PMAX', 'PMIN', 'QMAX', 'QMIN'}} ...
0998 );
0999 end
1000
1001
1002
1003 function slims = softlims_defaults(mpc, mpopt)
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
1015 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
1016 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
1017 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
1018 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
1019 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
1020 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
1021 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
1022
1023
1024 lims = softlims_lim2mat(mpopt);
1025 if isfield(mpc, 'softlims')
1026 slims = mpc.softlims;
1027 else
1028 slims = struct();
1029 end
1030 if isempty(mpopt)
1031 warning('softlims_defaults: Assuming ''mpopt.opf.softlims.default'' = 1, since mpopt was not provided.');
1032 use_default = 1;
1033 else
1034 use_default = mpopt.opf.softlims.default;
1035 end
1036
1037
1038 max_gen_cost = max(margcost(mpc.gencost, mpc.gen(:, PMAX)));
1039
1040
1041 for lm = fieldnames(lims).'
1042 lim = lm{:};
1043 mat = mpc.order.ext.(lims.(lim).mat_name);
1044 col = lims.(lim).col;
1045 specified = isfield(slims, lim);
1046 if ~specified && ~use_default
1047 slims.(lim).hl_mod = 'none';
1048 else
1049 if specified
1050 s = slims.(lim);
1051
1052
1053 if isfield(s, 'hl_mod') && strcmp(s.hl_mod, 'none')
1054 continue;
1055 end
1056 else
1057 s = struct();
1058 end
1059
1060
1061
1062
1063
1064 default_base_cost = max(1000, 2 * max_gen_cost);
1065 cost_scale = struct( ...
1066 'VMAX', 100, 'VMIN', 100, ...
1067 'ANGMAX', 1, 'ANGMIN', 1, ...
1068 'RATE_A', 1, ...
1069 'PMIN', 1, 'PMAX', 1, ...
1070 'QMIN', 1, 'QMAX', 1 );
1071 default_cost = cost_scale.(lim) * default_base_cost;
1072
1073
1074
1075 idxfull = softlims_default_idx(mat, lim, col);
1076
1077
1078 if ismember(lim, {'VMAX', 'VMIN'}) && ...
1079 (~isfield(s, 'idx') || isempty(s.idx)) && ...
1080 isfield(s, 'busnum') && ~isempty(s.busnum)
1081 s.idx = find(ismember(mat(:, BUS_I), s.busnum));
1082 end
1083
1084
1085 if isfield(s, 'idx')
1086 if size(s.idx, 2) > 1
1087 s.idx = s.idx.';
1088 if size(s.idx, 2) > 1
1089 error('softlim_defaults: mpc.softlims.%s.idx must be a scalar or vector', lim);
1090 end
1091 end
1092 else
1093 s.idx = idxfull;
1094 end
1095
1096
1097 if isempty(s.idx)
1098 s.hl_mod = 'none';
1099 slims.(lim) = s;
1100 continue;
1101 end
1102
1103
1104 if isfield(s, 'cost') && ~isempty(s.cost)
1105 if isscalar(s.cost)
1106 s.cost = s.cost * ones(size(s.idx));
1107 end
1108 else
1109 s.cost = default_cost * ones(size(s.idx));
1110 end
1111
1112
1113 ubsign = struct('VMAX', 1 , 'VMIN', -1, 'RATE_A', 1, ...
1114 'PMAX', 1, 'PMIN', -1, 'QMAX', 1, 'QMIN', -1, ...
1115 'ANGMAX', 1, 'ANGMIN', -1);
1116 shift_defaults = struct('VMAX', 0.25 , 'VMIN', 0.25, 'RATE_A', 10, ...
1117 'PMAX', 10, 'PMIN', 10, 'QMAX', 10, 'QMIN', 10, ...
1118 'ANGMAX', 10, 'ANGMIN', 10);
1119 if isfield(s, 'hl_mod')
1120 switch s.hl_mod
1121 case 'remove'
1122 case 'scale'
1123 if ~isfield(s, 'hl_val')
1124
1125
1126 orig_lim = mat(s.idx, col);
1127 s.hl_val = 2 * ones(size(s.idx));
1128 k = find(ubsign.(lim) * orig_lim < 0);
1129 s.hl_val(k) = 1/2;
1130 end
1131 case 'shift'
1132 if ~isfield(s, 'hl_val')
1133 s.hl_val = shift_defaults.(lim);
1134 end
1135 case 'replace'
1136 otherwise
1137 error('softlims_element_defaults: unknown hard limit modification %s', s.hl_mod);
1138 end
1139 else
1140 switch lim
1141 case 'PMIN'
1142
1143
1144 s.hl_mod = 'replace';
1145 s.hl_val = zeros(size(s.idx));
1146 s.hl_val(mat(s.idx, PMIN) < 0) = -Inf;
1147 case 'VMIN'
1148
1149 s.hl_mod = 'replace';
1150 s.hl_val = 0;
1151 case {'QMIN', 'ANGMIN'}
1152 s.hl_mod = 'remove';
1153 s.hl_val = -Inf;
1154 otherwise
1155 s.hl_mod = 'remove';
1156 s.hl_val = Inf;
1157 end
1158 end
1159
1160 slims.(lim) = s;
1161 end
1162 end
1163
1164
1165
1166 function idx = softlims_default_idx(mat, lim, col)
1167
1168
1169
1170
1171
1172
1173
1174 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
1175 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
1176 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
1177 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
1178 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
1179 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
1180 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
1181 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
1182
1183
1184
1185 switch lim
1186 case {'VMAX', 'VMIN'}
1187
1188 idx = [1:size(mat, 1)]';
1189 case {'ANGMAX', 'ANGMIN'}
1190
1191 idx = find(mat(:, BR_STATUS) > 0 & mat(:, col) & abs(mat(:, col)) < 360 );
1192 case 'RATE_A'
1193
1194 idx = find(mat(:, BR_STATUS) > 0 & mat(:, RATE_A) > 0);
1195 case {'PMAX', 'QMAX', 'QMIN'}
1196
1197 idx = find( (mat(:, GEN_STATUS) > 0) & ~isload(mat) & ~isinf(mat(:, col)) );
1198 case 'PMIN'
1199
1200 idx = find( (mat(:, GEN_STATUS) > 0) & ~isload(mat) );
1201 end
1202
1203
1204
1205 function slims = softlims_init(mpc, mpopt)
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
1217 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
1218 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
1219 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
1220 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
1221 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
1222 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
1223 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
1224
1225
1226 lims = softlims_lim2mat(mpopt);
1227 slims = mpc.softlims;
1228
1229
1230 for lm = fieldnames(lims).'
1231 lim = lm{:};
1232 s = slims.(lim);
1233
1234
1235 if isfield(s, 'hl_mod') && strcmp(s.hl_mod, 'none')
1236 continue;
1237 end
1238
1239 mat = mpc.order.ext.(lims.(lim).mat_name);
1240 col = lims.(lim).col;
1241 d = lims.(lim).dir;
1242
1243
1244
1245 idxfull = softlims_default_idx(mat, lim, col);
1246
1247
1248
1249 idxmask = ismember(s.idx, idxfull);
1250 s.idx = s.idx(idxmask);
1251
1252
1253 if isempty(s.idx)
1254 slims.(lim) = s;
1255 continue;
1256 end
1257
1258
1259 s.sav = mat(s.idx, col);
1260
1261
1262
1263 try
1264 s.cost = s.cost(idxmask);
1265 catch ME
1266 warning('softlims_init: something went wrong when handling the cost for limit %s. Perhaps the size of the ''cost'' vector didn''t match the size of the ''idx'' vector?', lim);
1267 rethrow(ME);
1268 end
1269
1270
1271
1272
1273 if isfield(s, 'hl_mod')
1274 switch s.hl_mod
1275 case 'remove'
1276 s.ub = Inf(size(s.idx));
1277 case 'scale'
1278
1279
1280 switch vectorcheck(s, idxmask)
1281 case 0
1282 error('softlims_init: provided hl_val vector does not conform to idx vector. When specifying hl_val as a vector idx must also be explicitly given');
1283 case 1
1284 s.ub = d * s.sav .* (s.hl_val(idxmask) - 1);
1285 case 2
1286 s.ub = d * s.sav * (s.hl_val - 1);
1287 case {1, 3}
1288 s.ub = d * s.sav .* (s.hl_val - 1);
1289 end
1290 case 'shift'
1291
1292
1293 switch vectorcheck(s, idxmask)
1294 case 0
1295 error('softlims_init: provided hl_val vector does not conform to idx vector. When specifying hl_val as a vector idx must also be explicitly given');
1296 case 1
1297 s.ub = s.hl_val(idxmask);
1298 case {2, 3}
1299 s.ub = s.hl_val * ones(size(s.idx));
1300 end
1301 case 'replace'
1302
1303
1304 switch vectorcheck(s, idxmask)
1305 case 0
1306 error('softlims_init: provided hl_val vector does not conform to idx vector. When specifying hl_val as a vector idx must also be explicitly given');
1307 case 1
1308 s.ub = d * ( s.hl_val(idxmask) - s.sav );
1309 case 2
1310 s.ub = d * ( s.hl_val - s.sav );
1311 case 3
1312 error('softlims_init: for hard limit ''replace'' modification, replacement value hl_val must be specified');
1313 end
1314 end
1315
1316 if any(s.ub < 0)
1317 error('softlims_init: some hl_val for %s results in more restrictive, rather than relaxed, hard constraint.', lim);
1318 end
1319 end
1320
1321
1322 switch lim
1323 case {'ANGMAX', 'ANGMIN', 'RATE_A'}
1324 s.rval = 0;
1325 case {'VMAX', 'PMAX', 'QMAX'}
1326 s.rval = Inf;
1327 case {'QMIN', 'PMIN','VMIN'}
1328 s.rval = -Inf;
1329 otherwise
1330 error('softlims_init: limit %s does not have an ''rval'' assigned', lim);
1331 end
1332 slims.(lim) = s;
1333 end
1334
1335
1336
1337 function out = vectorcheck(s, idx)
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347 if ~isfield(s, 'hl_val')
1348 out = 3;
1349 elseif isscalar(s.hl_val)
1350 out = 2;
1351 else
1352 out = all(size(s.hl_val) == size(idx));
1353 end
1354
1355
1356
1357 function [h, dh] = softlims_fcn(x, mpc, Yf, Yt, il, mpopt, Fmax)
1358
1359
1360
1361
1362 lim_type = upper(mpopt.opf.flow_lim(1));
1363
1364
1365
1366
1367 if nargout == 1
1368 h = opf_branch_flow_fcn(x(1:2), mpc, Yf, Yt, il, mpopt);
1369 else
1370 [h, dh] = opf_branch_flow_fcn(x(1:2), mpc, Yf, Yt, il, mpopt);
1371 end
1372 flv = x{3};
1373 if lim_type == 'P'
1374
1375
1376 tmp2 = Fmax + flv;
1377 else
1378
1379
1380 tmp1 = Fmax + flv;
1381 tmp2 = tmp1.^2;
1382 end
1383 h = h - [tmp2; tmp2];
1384 if nargout == 2
1385 ns = length(il);
1386 if lim_type == 'P'
1387 tmp3 = -speye(ns, ns);
1388 else
1389 tmp3 = spdiags(-2*tmp1, 0, ns, ns);
1390 end
1391 dh = [ dh [tmp3; tmp3] ];
1392 end
1393
1394
1395
1396 function d2H = softlims_hess(x, lambda, mpc, Yf, Yt, il, mpopt)
1397
1398
1399
1400
1401 lim_type = upper(mpopt.opf.flow_lim(1));
1402
1403
1404 d2H = opf_branch_flow_hess(x(1:2), lambda, mpc, Yf, Yt, il, mpopt);
1405 nh = size(d2H, 1);
1406 ns = length(il);
1407
1408 if lim_type == 'P'
1409 tmp = sparse(ns, ns);
1410 else
1411 tmp = spdiags(-2*lambda, 0, ns, ns);
1412 end
1413 d2H = [ d2H sparse(nh, ns);
1414 sparse(ns, nh) tmp ];