Skip to content

Commit

Permalink
Fix #45 where most does not properly handle cases with contingencies …
Browse files Browse the repository at this point in the history 
…defined only in some periods/scenarios.

- Check all periods/scenarios for contab before turning off contingencies when most.security_constrainted = -1
- Remove requirement for mdi.cont(1,1).contab to be non-empty when most.security_constrained = 1
- Update most_summary() to skip display of non-existent contingencies.
  • Loading branch information
@rdzman
rdzman committed May 22, 2024
1 parent 9a7d9d4 commit 7f616a4
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 49 deletions.
78 changes: 46 additions & 32 deletions lib/most.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,9 +144,19 @@
mo.IncludeFixedReserves = 0;
end
end
if mo.SecurityConstrained % check that at least some contingency is specified
have_contingency = 0;
for t = 1:nt
for j = 1:mdi.idx.nj(t)
if isfield(mdi, 'cont') && isfield(mdi.cont(t,j), 'contab') && ...
~isempty(mdi.cont(t,j).contab)
have_contingency = 1; % found a contingency
end
end
end
end
if mo.SecurityConstrained == -1
if isfield(mdi, 'cont') && isfield(mdi.cont(1,1), 'contab') && ...
~isempty(mdi.cont(1,1).contab)
if have_contingency
mo.SecurityConstrained = 1;
else
mo.SecurityConstrained = 0;
Expand All @@ -165,9 +175,8 @@
isfield(mdi.FixedReserves(1,1,1), 'req'))
error('most: MDI.FixedReserves(t,j,k) must be specified when MPOPT.most.fixed_res = 1');
end
if mo.SecurityConstrained && ~(isfield(mdi, 'cont') && ...
isfield(mdi.cont(1,1), 'contab') && ~isempty(mdi.cont(1,1).contab))
error('most: MDI.cont(t,j).contab cannot be empty when MPOPT.most.security_constraints = 1');
if mo.SecurityConstrained && ~have_contingency
error('most: MDI.cont(t,j).contab cannot be empty for all t, j when MPOPT.most.security_constraints = 1');
end
if mo.IncludeFixedReserves && mo.SecurityConstrained
warning('most: Using MPOPT.most.fixed_res = 1 and MPOPT.most.security_constraints = 1 together is not recommended.');
Expand Down Expand Up @@ -446,35 +455,40 @@
mdi.StepProb(t) = sum(scenario_probs); % probability of making it to the t-th step
if mdi.SecurityConstrained
for j = 1:mdi.idx.nj(t)
[tmp, ii] = sort(mdi.cont(t,j).contab(:, CT_LABEL)); %sort in ascending contingency label
contab = mdi.cont(t,j).contab(ii, :);
rowdecomlist = ones(size(contab,1), 1);
for l = 1:size(contab, 1)
if contab(l, CT_TABLE) == CT_TGEN && contab(l, CT_COL) == GEN_STATUS ...
&& contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % gen turned off
&& mdi.flow(t,j,1).mpc.gen(contab(l, CT_ROW), GEN_STATUS) <= 0 % but it was off on input
rowdecomlist(l) = 0;
elseif contab(l, CT_TABLE) == CT_TBRCH && contab(l, CT_COL) == BR_STATUS ...
&& contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % branch taken out
&& mdi.flow(t,j,1).mpc.branch(contab(l, CT_ROW), BR_STATUS) <= 0 % but it was off on input
rowdecomlist(l) = 0;
if isempty(mdi.cont(t,j).contab)
mdi.idx.nc(t, j) = 0;
mdi.CostWeights(1, j, t) = 1;
else
[tmp, ii] = sort(mdi.cont(t,j).contab(:, CT_LABEL)); %sort in ascending contingency label
contab = mdi.cont(t,j).contab(ii, :);
rowdecomlist = ones(size(contab,1), 1);
for l = 1:size(contab, 1)
if contab(l, CT_TABLE) == CT_TGEN && contab(l, CT_COL) == GEN_STATUS ...
&& contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % gen turned off
&& mdi.flow(t,j,1).mpc.gen(contab(l, CT_ROW), GEN_STATUS) <= 0 % but it was off on input
rowdecomlist(l) = 0;
elseif contab(l, CT_TABLE) == CT_TBRCH && contab(l, CT_COL) == BR_STATUS ...
&& contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % branch taken out
&& mdi.flow(t,j,1).mpc.branch(contab(l, CT_ROW), BR_STATUS) <= 0 % but it was off on input
rowdecomlist(l) = 0;
end
end
contab = contab(rowdecomlist ~= 0, :);
mdi.cont(t, j).contab = contab;
clist = unique(contab(:, CT_LABEL));
mdi.idx.nc(t, j) = length(clist);
k = 2;
for label = clist'
mdi.flow(t, j, k).mpc = apply_changes(label, mdi.flow(t, j, 1).mpc, contab);
ii = find( label == contab(:, CT_LABEL) );
mdi.CostWeights(k, j, t) = contab(ii(1), CT_PROB);
mdi.idx.nb(t, j, k) = size(mdi.flow(t, j, k).mpc.bus, 1);
mdi.idx.ny(t, j, k) = length(find(mdi.flow(t, j, 1).mpc.gencost(:, MODEL) == PW_LINEAR));
k = k + 1;
end
mdi.CostWeights(1, j, t) = 1 - sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1, j, t));
mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t) = scenario_probs(j) * mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t);
end
contab = contab(rowdecomlist ~= 0, :);
mdi.cont(t, j).contab = contab;
clist = unique(contab(:, CT_LABEL));
mdi.idx.nc(t, j) = length(clist);
k = 2;
for label = clist'
mdi.flow(t, j, k).mpc = apply_changes(label, mdi.flow(t, j, 1).mpc, contab);
ii = find( label == contab(:, CT_LABEL) );
mdi.CostWeights(k, j, t) = contab(ii(1), CT_PROB);
mdi.idx.nb(t, j, k) = size(mdi.flow(t, j, k).mpc.bus, 1);
mdi.idx.ny(t, j, k) = length(find(mdi.flow(t, j, 1).mpc.gencost(:, MODEL) == PW_LINEAR));
k = k + 1;
end
mdi.CostWeights(1, j, t) = 1 - sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1, j, t));
mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t) = scenario_probs(j) * mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t);
end
else
for j = 1:mdi.idx.nj(t)
Expand Down
37 changes: 20 additions & 17 deletions lib/most_summary.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,28 +68,30 @@
nl = size(mpc.branch, 1);
ng = size(mpc.gen, 1);
nt = mdo.idx.nt;
nj_max = max(mdo.idx.nj);
nc_max = max(max(mdo.idx.nc));
nj = mdo.idx.nj;
nc = mdo.idx.nc;
nj_max = max(nj);
nc_max = max(max(nc));
ns = mdo.idx.ns;

%% summarize results
psi = zeros(nt, nj_max, nc_max+1);
Pg = zeros(ng, nt, nj_max, nc_max+1);
Pd = zeros(nb, nt, nj_max, nc_max+1);
psi = NaN(nt, nj_max, nc_max+1);
Pg = NaN(ng, nt, nj_max, nc_max+1);
Pd = NaN(nb, nt, nj_max, nc_max+1);
if mdo.idx.ntramp
Rup = mdo.results.Rrp;
Rdn = mdo.results.Rrm;
else
Rup = [];
Rdn = [];
end
u = zeros(ng, nt);
lamP = zeros(nb, nt, nj_max, nc_max+1);
muF = zeros(nl, nt, nj_max, nc_max+1);
Pf = zeros(nl, nt, nj_max, nc_max+1);
u = NaN(ng, nt);
lamP = NaN(nb, nt, nj_max, nc_max+1);
muF = NaN(nl, nt, nj_max, nc_max+1);
Pf = NaN(nl, nt, nj_max, nc_max+1);
for t = 1:nt
for j = 1:mdo.idx.nj(t)
for k = 1:mdo.idx.nc(t,j)+1
for j = 1:nj(t)
for k = 1:nc(t,j)+1
rr = mdo.flow(t,j,k).mpc;
psi(t, j, k) = mdo.CostWeightsAdj(k, j, t);
u(:, t) = rr.gen(:, GEN_STATUS);
Expand Down Expand Up @@ -159,19 +161,19 @@
end
fprintf('\n');

print_most_summary_section('PG', 'Gen', nt, nj_max, nc_max, Pg);
print_most_summary_section('PG', 'Gen', nt, nj_max, nc, Pg);
if mdo.idx.ntramp
print_most_summary_section('RAMP UP', 'Gen', nt, 1, 0, Rup);
print_most_summary_section('RAMP DOWN', 'Gen', nt, 1, 0, Rdn);
end
print_most_summary_section('FIXED LOAD', 'Bus', nt, nj_max, nc_max, Pd);
print_most_summary_section('FIXED LOAD', 'Bus', nt, nj_max, nc, Pd);
if ns
print_most_summary_section('ESS E[SoC]', 'ESS', nt, 1, 0, SoC);
end
if mdo.DCMODEL
print_most_summary_section('LAM_P', 'Bus', nt, nj_max, nc_max, lamP);
print_most_summary_section('PF', 'Brch', nt, nj_max, nc_max, Pf);
print_most_summary_section('MU_F', 'Brch', nt, nj_max, nc_max, muF);
print_most_summary_section('LAM_P', 'Bus', nt, nj_max, nc, lamP);
print_most_summary_section('PF', 'Brch', nt, nj_max, nc, Pf);
print_most_summary_section('MU_F', 'Brch', nt, nj_max, nc, muF);
end
end

Expand All @@ -180,7 +182,7 @@
end

%%---------------------------------------------------------
function print_most_summary_section(label, section_type, nt, nj_max, nc_max, data, tol)
function print_most_summary_section(label, section_type, nt, nj_max, nc, data, tol)
if nargin < 7
tol = 1e-4;
end
Expand All @@ -189,6 +191,7 @@ function print_most_summary_section(label, section_type, nt, nj_max, nc_max, dat
fprintf('\n==========%-12s==========\n', sprintf('%s%s', bl, label));
if any(data(:))
for j = 1:nj_max
nc_max = max(nc(:,j));
for k = 1:nc_max+1
if nj_max > 1 || nc_max > 0
fprintf('\nSCENARIO %d', j);
Expand Down

0 comments on commit 7f616a4

Please sign in to comment.