Skip to content

Commit

Permalink
Use opt_model/get_soln() to extract results
Browse files Browse the repository at this point in the history 
rather than doing the indexing manually.
  • Loading branch information
@rdzman
rdzman committed Sep 9, 2020
1 parent a43f587 commit 1cf19d4
Show file tree
Hide file tree
Showing 3 changed files with 54 additions and 49 deletions.
9 changes: 7 additions & 2 deletions CHANGES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,15 @@ Change history for MOST
Changes since 1.0.2
-------------------

*Requires MATPOWER with 7.1 or later (for MP-Opt-Model 2.2 or later).*

#### 9/9/20
- Use `@opt_model/get_soln()` to extract variable and shadow price
results, rather than doing the indexing manually.

#### 3/19/20
- Convert to using `@opt_model/solve()` method rather than calling
`miqps_matpower()` or `qps_matpower()` directly. Requires MATPOWER
7.1 or later.
`miqps_matpower()` or `qps_matpower()` directly.
- **INCOMPATIBLE CHANGE**: Update objective function value returned in
`mdo.QP.f` to include the previously missing constant term.

Expand Down
1 change: 1 addition & 0 deletions docs/src/MOST-manual/MOST-manual.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3353,6 +3353,7 @@ \subsubsection*{Other Changes}
\item Update \code{most\_summary()} to include sections for fixed loads
and storage expected stored energy.
\item Convert to using \code{@opt\_model/solve()} method rather than calling \code{miqps\_matpower()} or \code{qps\_matpower()} directly. Requires \matpower{} \hl{7.1} or later.
\item Use \code{@opt\_model/get\_soln()} to extract variable and shadow price results, rather than doing the indexing manually. Requires \matpower{} \hl{7.1} or later.
\end{itemize}

\subsubsection*{Incompatible Changes}
Expand Down
93 changes: 46 additions & 47 deletions lib/most.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1980,7 +1980,6 @@
mdo = mdi; %% initialize output
mdo.om = om.copy(); %% make copy of opt_model object, so changes to
%% output obj (mdo) don't modify input obj (mdi)
[vv, ll] = mdo.om.get_idx();
if mpopt.most.solve_model
%% check consistency of model options (in case mdi was built in previous call)
if mdi.DCMODEL ~= mo.DCMODEL
Expand Down Expand Up @@ -2037,9 +2036,10 @@
fprintf('- Post-processing results.\n');
end
if success
om = mdo.om;
for t = 1:nt
if UC
mdo.UC.CommitSched(:, t) = mdo.QP.x(vv.i1.u(t):vv.iN.u(t));
mdo.UC.CommitSched(:, t) = om.get_soln('var', 'u', {t});
end
for j = 1:mdi.idx.nj(t)
for k = 1:mdi.idx.nc(t,j)+1
Expand All @@ -2049,7 +2049,7 @@
mpc.bus(:, VM) = 1;
end
% Injections and shadow prices
mpc.gen(:, PG) = baseMVA * mdo.QP.x(vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k));
mpc.gen(:, PG) = baseMVA * om.get_soln('var', 'Pg', {t,j,k});
%% need to update Qg for loads consistent w/constant power factor
Pmin = mpc.gen(:, PMIN);
Qmin = mpc.gen(:, QMIN);
Expand All @@ -2059,10 +2059,11 @@
mpc.gen(ivl, QG) = mpc.gen(ivl, PG) .* Qlim ./ Pmin(ivl);
if mdo.DCMODEL
%% bus angles
mpc.bus(:, VA) = (180/pi) * mdo.QP.x(vv.i1.Va(t,j,k):vv.iN.Va(t,j,k));
Va = om.get_soln('var', 'Va', {t,j,k});
mpc.bus(:, VA) = (180/pi) * Va;

%% nodal prices
price = (mdo.QP.lambda.mu_u(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))-mdo.QP.lambda.mu_l(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))) / baseMVA;
price = (om.get_soln('lin', 'mu_u', 'Pmis', {t,j,k})-om.get_soln('lin', 'mu_l', 'Pmis', {t,j,k})) / baseMVA;
mpc.bus(:, LAM_P) = price;

%% line flows and line limit shadow prices
Expand All @@ -2073,50 +2074,48 @@
mpc.branch(:, MU_SF) = 0;
mpc.branch(:, MU_ST) = 0;
ion = find(mpc.branch(:, BR_STATUS));
rows = ll.i1.Pf(t,j,k):ll.iN.Pf(t,j,k);
cols = vv.i1.Va(t,j,k):vv.iN.Va(t,j,k);
lf = baseMVA * (mdo.QP.A(rows,cols) * mdo.QP.x(cols) + mdo.flow(t,j,k).PLsh);
AA = om.params_lin_constraint('Pf', {t,j,k});
lf = baseMVA * (AA * Va + mdo.flow(t,j,k).PLsh);
mpc.branch(ion, PF) = lf;
mpc.branch(ion, PT) = -lf;
mpc.branch(ion, MU_SF) = mdo.QP.lambda.mu_u(rows) / baseMVA;
mpc.branch(ion, MU_ST) = mdo.QP.lambda.mu_l(rows) / baseMVA;
mpc.branch(ion, MU_SF) = om.get_soln('lin', 'mu_u', 'Pf', {t,j,k}) / baseMVA;
mpc.branch(ion, MU_ST) = om.get_soln('lin', 'mu_l', 'Pf', {t,j,k}) / baseMVA;
else
%% system price
price = (mdo.QP.lambda.mu_l(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))-mdo.QP.lambda.mu_u(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))) / baseMVA;
price = (om.get_soln('lin', 'mu_l', 'Pmis', {t,j,k})-om.get_soln('lin', 'mu_u', 'Pmis', {t,j,k})) / baseMVA;
mpc.bus(:, LAM_P) = price;
end
if UC
% igenon does not contain gens ousted because of a contingency or
% a forced-off UC.CommitKey
igenon = find(mpc.gen(:, GEN_STATUS));
u = mdo.QP.x(vv.i1.u(t):vv.iN.u(t));
mpc.gen(igenon, GEN_STATUS) = u(igenon);
mpc.gen(igenon, GEN_STATUS) = mdo.UC.CommitSched(igenon, t);
gs = mpc.gen(igenon, GEN_STATUS) > 0; % gen status
mpc.gen(:, MU_PMAX) = 0;
mpc.gen(:, MU_PMIN) = 0;
mpc.gen(igenon, MU_PMAX) = gs .* ...
mdo.QP.lambda.mu_u(ll.i1.uPmax(t,j,k):ll.iN.uPmax(t,j,k)) / baseMVA;
om.get_soln('lin', 'mu_u', 'uPmax', {t,j,k}) / baseMVA;
mpc.gen(igenon, MU_PMIN) = gs .* ...
mdo.QP.lambda.mu_u(ll.i1.uPmin(t,j,k):ll.iN.uPmin(t,j,k)) / baseMVA;
om.get_soln('lin', 'mu_u', 'uPmin', {t,j,k}) / baseMVA;
if mdo.QCoordination
mpc.gen(:, MU_QMAX) = 0;
mpc.gen(:, MU_QMIN) = 0;
mpc.gen(igenon, MU_QMAX) = gs .* ...
mdo.QP.lambda.mu_u(ll.i1.uQmax(t,j,k):ll.iN.uQmax(t,j,k)) / baseMVA;
om.get_soln('lin', 'mu_u', 'uQmax', {t,j,k}) / baseMVA;
mpc.gen(igenon, MU_QMIN) = gs .* ...
mdo.QP.lambda.mu_u(ll.i1.uQmin(t,j,k):ll.iN.uQmin(t,j,k)) / baseMVA;
om.get_soln('lin', 'mu_u', 'uQmin', {t,j,k}) / baseMVA;
end
else
gs = mpc.gen(:, GEN_STATUS) > 0; % gen status
mpc.gen(:, MU_PMAX) = gs .* ...
mdo.QP.lambda.upper(vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k)) / baseMVA;
om.get_soln('var', 'mu_u', 'Pg', {t,j,k}) / baseMVA;
mpc.gen(:, MU_PMIN) = gs .* ...
mdo.QP.lambda.lower(vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k)) / baseMVA;
om.get_soln('var', 'mu_l', 'Pg', {t,j,k}) / baseMVA;
if mdo.QCoordination
mpc.gen(:, MU_QMAX) = gs .* ...
mdo.QP.lambda.upper(vv.i1.Qg(t,j,k):vv.iN.Qg(t,j,k)) / baseMVA;
om.get_soln('var', 'mu_u', 'Qg', {t,j,k}) / baseMVA;
mpc.gen(:, MU_QMIN) = gs .* ...
mdo.QP.lambda.lower(vv.i1.Qg(t,j,k):vv.iN.Qg(t,j,k)) / baseMVA;
om.get_soln('var', 'mu_l', 'Qg', {t,j,k}) / baseMVA;
end
end
if mdi.IncludeFixedReserves
Expand All @@ -2125,34 +2124,34 @@
r.R = z;
r.prc = z;
r.mu = struct('l', z, 'u', z, 'Pmax', z);
r.totalcost = sum(mdo.om.eval_quad_cost(mdo.QP.x, 'Rcost', {t,j,k}));
r.R(r.igr) = mdo.QP.x(vv.i1.R(t,j,k):vv.iN.R(t,j,k)) * baseMVA;
r.totalcost = sum(om.get_soln('qdc', 'Rcost', {t,j,k}));
r.R(r.igr) = om.get_soln('var', 'R', {t,j,k}) * baseMVA;
R_mu_l = om.get_soln('lin', 'mu_l', 'Rreq', {t,j,k});
for gg = r.igr
iz = find(r.zones(:, gg));
kk = ll.i1.Rreq(t,j,k):ll.iN.Rreq(t,j,k);
r.prc(gg) = sum(mdo.QP.lambda.mu_l(kk(iz))) / baseMVA;
r.prc(gg) = sum(R_mu_l(iz)) / baseMVA;
end
r.mu.l(r.igr) = mdo.QP.lambda.lower(vv.i1.R(t,j,k):vv.iN.R(t,j,k)) / baseMVA;
r.mu.u(r.igr) = mdo.QP.lambda.upper(vv.i1.R(t,j,k):vv.iN.R(t,j,k)) / baseMVA;
r.mu.Pmax(r.igr) = mdo.QP.lambda.mu_u(ll.i1.Pg_plus_R(t,j,k):ll.iN.Pg_plus_R(t,j,k)) / baseMVA;
r.mu.l(r.igr) = om.get_soln('var', 'mu_l', 'R', {t,j,k}) / baseMVA;
r.mu.u(r.igr) = om.get_soln('var', 'mu_u', 'R', {t,j,k}) / baseMVA;
r.mu.Pmax(r.igr) = om.get_soln('lin', 'mu_u', 'Pg_plus_R', {t,j,k}) / baseMVA;
mpc.reserves = r;
end
mdo.flow(t,j,k).mpc = mpc; %% stash modified mpc in output struct
end
end
% Contract, contingency reserves, energy limits
mdo.results.Pc(:,t) = baseMVA * mdo.QP.x(vv.i1.Pc(t):vv.iN.Pc(t));
mdo.results.Rpp(:,t) = baseMVA * mdo.QP.x(vv.i1.Rpp(t):vv.iN.Rpp(t));
mdo.results.Rpm(:,t) = baseMVA * mdo.QP.x(vv.i1.Rpm(t):vv.iN.Rpm(t));
mdo.results.Pc(:,t) = baseMVA * om.get_soln('var', 'Pc', {t});
mdo.results.Rpp(:,t) = baseMVA * om.get_soln('var', 'Rpp', {t});
mdo.results.Rpm(:,t) = baseMVA * om.get_soln('var', 'Rpm', {t});
if ns
mdo.results.Sm(:,t) = baseMVA * mdo.QP.x(vv.i1.Sm(t):vv.iN.Sm(t));
mdo.results.Sp(:,t) = baseMVA * mdo.QP.x(vv.i1.Sp(t):vv.iN.Sp(t));
mdo.results.Sm(:,t) = baseMVA * om.get_soln('var', 'Sm', {t});
mdo.results.Sp(:,t) = baseMVA * om.get_soln('var', 'Sp', {t});
end
end
% Ramping reserves
for t = 1:mdo.idx.ntramp
mdo.results.Rrp(:,t) = baseMVA * mdo.QP.x(vv.i1.Rrp(t):vv.iN.Rrp(t));
mdo.results.Rrm(:,t) = baseMVA * mdo.QP.x(vv.i1.Rrm(t):vv.iN.Rrm(t));
mdo.results.Rrp(:,t) = baseMVA * om.get_soln('var', 'Rrp', {t});
mdo.results.Rrm(:,t) = baseMVA * om.get_soln('var', 'Rrm', {t});
end
% Expected energy prices for generators, per generator and per period,
% both absolute and conditional on making it to that period
Expand All @@ -2172,13 +2171,13 @@
mdo.results.RppPrices = zeros(ng, nt);
mdo.results.RpmPrices = zeros(ng, nt);
for t = 1:nt
mdo.results.RppPrices(:, t) = mdo.QP.lambda.lower(vv.i1.Rpp(t):vv.iN.Rpp(t)) / baseMVA;
mdo.results.RpmPrices(:, t) = mdo.QP.lambda.lower(vv.i1.Rpm(t):vv.iN.Rpm(t)) / baseMVA;
mdo.results.RppPrices(:, t) = om.get_soln('var', 'mu_l', 'Rpp', {t}) / baseMVA;
mdo.results.RpmPrices(:, t) = om.get_soln('var', 'mu_l', 'Rpm', {t}) / baseMVA;
for j = 1:mdi.idx.nj(t);
for k = 1:mdi.idx.nc(t,j)+1
ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
mdo.results.RppPrices(ii, t) = mdo.results.RppPrices(ii, t) + mdo.QP.lambda.mu_l(ll.i1.dPpRp(t,j,k):ll.iN.dPpRp(t,j,k)) / baseMVA;
mdo.results.RpmPrices(ii, t) = mdo.results.RpmPrices(ii, t) + mdo.QP.lambda.mu_l(ll.i1.dPmRm(t,j,k):ll.iN.dPmRm(t,j,k)) / baseMVA;
mdo.results.RppPrices(ii, t) = mdo.results.RppPrices(ii, t) + om.get_soln('lin', 'mu_l', 'dPpRp', {t,j,k}) / baseMVA;
mdo.results.RpmPrices(ii, t) = mdo.results.RpmPrices(ii, t) + om.get_soln('lin', 'mu_l', 'dPmRm', {t,j,k}) / baseMVA;
end
end
end
Expand All @@ -2190,17 +2189,17 @@
for j1 = 1:mdo.idx.nj(t)
for j2 = 1:mdo.idx.nj(t+1)
if mdi.tstep(t+1).TransMask(j2,j1)
mdo.results.RrpPrices(:, t) = mdo.results.RrpPrices(:, t) + mdo.QP.lambda.mu_l(ll.i1.Rrp(t,j1,j2):ll.iN.Rrp(t,j1,j2)) / baseMVA;
mdo.results.RrmPrices(:, t) = mdo.results.RrmPrices(:, t) + mdo.QP.lambda.mu_l(ll.i1.Rrm(t,j1,j2):ll.iN.Rrm(t,j1,j2)) / baseMVA;
mdo.results.RrpPrices(:, t) = mdo.results.RrpPrices(:, t) + om.get_soln('lin', 'mu_l', 'Rrp', {t,j1,j2}) / baseMVA;
mdo.results.RrmPrices(:, t) = mdo.results.RrmPrices(:, t) + om.get_soln('lin', 'mu_l', 'Rrm', {t,j1,j2}) / baseMVA;
end
end
end
end
% then last period only if specified for with terminal state
if ~mdo.OpenEnded
for j1 = 1:mdo.idx.nj(nt)
mdo.results.RrpPrices(:, nt) = mdo.results.RrpPrices(:, nt) + mdo.QP.lambda.mu_l(ll.i1.Rrp(nt,j1,1):ll.iN.Rrp(nt,j1,1)) / baseMVA;
mdo.results.RrmPrices(:, nt) = mdo.results.RrmPrices(:, nt) + mdo.QP.lambda.mu_l(ll.i1.Rrm(nt,j1,1):ll.iN.Rrm(nt,j1,1)) / baseMVA;
mdo.results.RrpPrices(:, nt) = mdo.results.RrpPrices(:, nt) + om.get_soln('lin', 'mu_l', 'Rrp', {nt,j1,1}) / baseMVA;
mdo.results.RrmPrices(:, nt) = mdo.results.RrmPrices(:, nt) + om.get_soln('lin', 'mu_l', 'Rrm', {nt,j1,1}) / baseMVA;
end
end
% Expected wear and tear costs per gen and period
Expand Down Expand Up @@ -2240,7 +2239,7 @@
end
% If Cyclic storage, pull InitialStorage value out of x
if ns && mdo.Storage.ForceCyclicStorage
mdo.Storage.InitialStorage = baseMVA * mdo.QP.x(vv.i1.S0:vv.iN.S0);
mdo.Storage.InitialStorage = baseMVA * om.get_soln('var', 'S0');
end
% Compute expected storage state trajectory
mdo.Storage.ExpectedStorageState = zeros(ns,nt);
Expand All @@ -2267,14 +2266,14 @@
if nzds
mdo.results.Z = zeros(nzds, ntds);
for t = 1:ntds
mdo.results.Z(:,t) = mdo.QP.x(vv.i1.Z(t):vv.iN.Z(t));
mdo.results.Z(:,t) = om.get_soln('var', 'Z', {t});
end
end
mdo.results.Y = zeros(nyds, ntds);
if nyds
for t = 1:ntds
mdo.results.Y(:, t) = ...
mdo.QP.A(ll.i1.DSy(t):ll.iN.DSy(t), :) * mdo.QP.x;
AA = om.params_lin_constraint('DSy', {t});
mdo.results.Y(:, t) = AA * mdo.QP.x;
end
end
end
Expand Down

0 comments on commit 1cf19d4

Please sign in to comment.