diff --git a/lib/most.m b/lib/most.m index 0c3cd99..b085fcc 100644 --- a/lib/most.m +++ b/lib/most.m @@ -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; @@ -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.'); @@ -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) diff --git a/lib/most_summary.m b/lib/most_summary.m index 0b62e09..804b26f 100644 --- a/lib/most_summary.m +++ b/lib/most_summary.m @@ -68,14 +68,16 @@ 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; @@ -83,13 +85,13 @@ 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); @@ -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 @@ -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 @@ -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);