From 8ee608bf109a598f08ed412b31743c1ec836ab07 Mon Sep 17 00:00:00 2001 From: Ray Zimmerman Date: Wed, 25 Oct 2023 12:49:33 -0600 Subject: [PATCH] Correct mistake in handling of ramp wear & tear for Detla_T ~= 1 --- CHANGES.md | 2 +- docs/src/MOST-manual/MOST-manual.tex | 5 ++-- lib/most.m | 42 ++++++++++++++-------------- lib/mostver.m | 2 +- 4 files changed, 26 insertions(+), 25 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 16f8d41..97df969 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -5,7 +5,7 @@ Change history for MOST since 1.2 --------- -#### 10/24/23 +#### 10/25/23 - Fix [issue #39][10] in which the value of `mdi.Delta_T`, the number of hours represented by each period, was not being accounted for in most of the terms in the objective function. diff --git a/docs/src/MOST-manual/MOST-manual.tex b/docs/src/MOST-manual/MOST-manual.tex index 28d6ff9..90c93ea 100644 --- a/docs/src/MOST-manual/MOST-manual.tex +++ b/docs/src/MOST-manual/MOST-manual.tex @@ -1469,7 +1469,7 @@ \subsubsection{{\tt xgd} -- Extra Generator Data ({\tt xGenData})} \item [*] {All fields are $n_g \times 1$ vectors of per-generator values.} \item [\dag] {These are defaults provided by \code{loadxgendata}. If \code{gen} is provided, either directly or as the \code{gen} field of \mpc{}, then \code{P = gen(:, PG)}, \code{C = gen(:, GEN\_STATUS)} and \code{R = 2*(gen(:, PMAX) - MIN(0, gen(:, PMIN)))}, otherwise $\code{C}~=~1$, $\code{R}~=~0$ and no default is provided for \code{P} (corresponding field is not optional).} \item [\dag\dag] {Ramping costs/restrictions from the initial dispatch at $t=0$ are ignored for single-period problems.} - \item [\ddag] {Each of these costs $C(\cdot)$ is presented in the formulation as a general function, but is implmented as a simple linear function of the form $C(x) = a x$, where the linear coefficient being supplied is $a$. The only exception is the ramping cost, which has the quadratic form $C_\delta^i(x) = a x^2$.} + \item [\ddag] {Each of these costs $C(\cdot)$ is presented in the formulation as a general function, but is implmented as a simple linear function of the form $C(x) = a x$, where the linear coefficient being supplied is $a$. The only exception is the ramping cost, which has the quadratic form $C_\delta^i(x) = \frac{1}{\Delta^2} a x^2$.} \item [\S] {Requires that \code{CommitKey} be present and non-empty.} \item [\P] {Sign is based on \code{C}\tnote{\dag}, i.e. $+\infty$ for \code{C} = 1, $-\infty$ for \code{C} = 0.} \end{tablenotes} @@ -1722,7 +1722,7 @@ \subsubsection{Input Data} \code{mpc} & I & & base system data, standard \matpower{} case struct\tnote{\ddag}, with \baseMVA{}, \bus{}, \gen{}, \branch{} and \gencost{} fields \\ \code{offer(t)} & I & & struct with offer data for period~$t$, see Table~\ref{tab:md_inputoffer} for details of sub-fields \\ \code{OpenEnded} & I & 1 & ignore terminal dispatch ramp constraints, \emph{deprecated} \\ -\code{RampWearCostCoeff(i,t)} & I & 0 & $n_g \times n_t$, cost of ramping of generator~$i$ from period~$t-1$ to $t$, coefficient $C_\delta^i$ for square of dispatch difference in \eqref{eq:rampcost} \\ +\code{RampWearCostCoeff(i,t)} & I & 0 & $n_g \times n_t$, cost of ramping of generator~$i$ from period~$t-1$ to $t$, coefficient $C_\delta^i$ for square of dispatch difference in \eqref{eq:rampcost}\tnote{\S} \\ \code{Storage} & B & & struct with parameters for storage units, see Table~\ref{tab:md_inputstorage} for the input fields \\ \code{TerminalPg(i)} & I & & $n_g \times 1$, injection of generator~$i$ at $t = n_t$, \emph{deprecated, untested} \\ \code{tstep(t)} & B & & $n_t \times 1$ struct of parameters related to period~$t$ \\ @@ -1736,6 +1736,7 @@ \subsubsection{Input Data} \item [*] {I = input, O = output, B = both, opt = taken from \matpower{} options.} \item [\dag] {See Section~\ref{sec:contab} for details. Note that, while \code{loadmd} assigns the same \code{contab} to all $t$ and $j$, it is possible to set different \code{contab} values manually and they will be respected by \code{most}.} \item [\ddag] {See Appendix~\ref{MUM-app:caseformat} in the \mum{} for details.} + \item [\S] {More precisely, $C_\delta^i(x) = \frac{1}{\Delta^2} a x^2$}, where $a$ is the corresponding value of \code{RampWearCostCoeff(i,t)}. \end{tablenotes} \end{threeparttable} \end{table} diff --git a/lib/most.m b/lib/most.m index edc77ec..c733079 100644 --- a/lib/most.m +++ b/lib/most.m @@ -1722,10 +1722,10 @@ om.init_indexed_name('qdc', 'RampWear', {nt+1, nj_max, nj_max}); end for j = 1:mdi.idx.nj(1) - w = mdi.tstep(1).TransMat(j,1); % the probability of going from initial state to jth - Q = mdi.Delta_T * spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,1), 0, ng, ng); - c = -mdi.Delta_T * w * baseMVA * mdi.RampWearCostCoeff(:,1) .* mdi.InitialPg; - k0 = mdi.Delta_T * w * 0.5 * mdi.RampWearCostCoeff(:,1)' * mdi.InitialPg.^2; + w = mdi.Delta_T * mdi.tstep(1).TransMat(j,1); % the probability of going from initial state to jth + Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,1) / mdi.Delta_T^2, 0, ng, ng); + c = -w * baseMVA * mdi.RampWearCostCoeff(:,1) .* mdi.InitialPg; + k0 = w * 0.5 * mdi.RampWearCostCoeff(:,1)' * mdi.InitialPg.^2; vs = struct('name', {'Pg'}, 'idx', {{1,j,1}}); om.add_quad_cost('RampWear', {1,j,1}, Q, c, k0, vs); end @@ -1733,8 +1733,8 @@ for t = 2:nt for j2 = 1:mdi.idx.nj(t) for j1 = 1:mdi.idx.nj(t-1) - w = mdi.tstep(t).TransMat(j2,j1) * mdi.CostWeights(1, j1, t-1); - h = mdi.Delta_T * w * baseMVA^2 * mdi.RampWearCostCoeff(:,t); + w = mdi.Delta_T * mdi.tstep(t).TransMat(j2,j1) * mdi.CostWeights(1, j1, t-1); + h = w * baseMVA^2 * mdi.RampWearCostCoeff(:,t) / mdi.Delta_T^2; i = (1:ng)'; j = ng+(1:ng)'; Q = sparse([i;j;i;j], [i;i;j;j], [h;-h;-h;h], 2*ng, 2*ng); @@ -1749,10 +1749,10 @@ % that makes sense for nt+1; all other fields in mdi.tstep(nt+1) can be empty. if ~mdi.OpenEnded for j = 1:mdi.idx.nj(nt) - w = mdi.tstep(nt+1).TransMat(1, j) * mdi.CostWeights(1, j, nt); - Q = mdi.Delta_T * spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,nt+1), 0, ng, ng); - c = -mdi.Delta_T * w * baseMVA * mdi.RampWearCostCoeff(:,nt+1) .* mdi.TerminalPg; - k0 = mdi.Delta_T * w * 0.5 * mdi.RampWearCostCoeff(:,nt+1)' * mdi.TerminalPg.^2; + w = mdi.Delta_T * mdi.tstep(nt+1).TransMat(1, j) * mdi.CostWeights(1, j, nt); + Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,nt+1) / mdi.Delta_T^2, 0, ng, ng); + c = -w * baseMVA * mdi.RampWearCostCoeff(:,nt+1) .* mdi.TerminalPg; + k0 = w * 0.5 * mdi.RampWearCostCoeff(:,nt+1)' * mdi.TerminalPg.^2; vs = struct('name', {'Pg'}, 'idx', {{nt,j,1}}); om.add_quad_cost('RampWear', {nt+1,j,1}, Q, c, k0, vs); end @@ -1772,7 +1772,7 @@ for t = 1:nt for j = 1:mdi.idx.nj(t) for k = 1:mdi.idx.nc(t,j)+1 - w = mdi.CostWeightsAdj(k,j,t); %% NOTE (k,j,t) order !!! + w = mdi.Delta_T * mdi.CostWeightsAdj(k,j,t); %% NOTE (k,j,t) order !!! % weighted polynomial energy costs for committed units gc = mdi.flow(t,j,k).mpc.gencost; @@ -1784,15 +1784,15 @@ if ncost > 3 error('most: polynomial generator costs of order higher than quadratic not supported'); elseif ncost == 3 - Q = sparse(ipol, ipol, mdi.Delta_T * 2 * w * baseMVA^2*gc(ipol, COST), ng, ng); + Q = sparse(ipol, ipol, 2 * w * baseMVA^2*gc(ipol, COST), ng, ng); else Q = sparse(ng,ng); end c = zeros(ng, 1); if ncost >= 2 - c(ipol) = mdi.Delta_T * w * baseMVA*gc(ipol, COST+ncost-2); + c(ipol) = w * baseMVA*gc(ipol, COST+ncost-2); end - k0 = mdi.Delta_T * w * sum(gc(ipol, COST+ncost-1)); + k0 = w * sum(gc(ipol, COST+ncost-1)); else %% non-uniform order of polynomials %% use a loop Q = sparse(ng,ng); @@ -1802,12 +1802,12 @@ if ncost > 3 error('most: polynomial generator costs of order higher than quadratic not supported'); elseif ncost == 3 - Q(i,i) = mdi.Delta_T * 2 * w * baseMVA^2*gc(i, COST); + Q(i,i) = 2 * w * baseMVA^2*gc(i, COST); end if ncost >= 2 - c(i) = mdi.Delta_T * w * baseMVA*gc(i, COST+ncost-2); + c(i) = w * baseMVA*gc(i, COST+ncost-2); end - k0 = mdi.Delta_T * w * gc(i, COST+ncost-1); + k0 = w * gc(i, COST+ncost-1); end end vs = struct('name', {'Pg'}, 'idx', {{t,j,k}}); @@ -1817,22 +1817,22 @@ % weighted y-variables for piecewise linear energy costs for committed units % ipwl = find( (mdi.flow(t,j,k).mpc.gen(:,GEN_STATUS) > 0) & (gc(:,MODEL) == PW_LINEAR)); if mdi.idx.ny(t,j,k) - c = mdi.Delta_T * w * ones(mdi.idx.ny(t,j,k),1); + c = w * ones(mdi.idx.ny(t,j,k),1); vs = struct('name', {'y'}, 'idx', {{t,j,k}}); om.add_quad_cost('Cy', {t,j,k}, [], c, 0, vs); end % inc and dec offers for each flow - c = mdi.Delta_T * w * baseMVA * mdi.offer(t).PositiveActiveDeltaPrice(:); + c = w * baseMVA * mdi.offer(t).PositiveActiveDeltaPrice(:); vs = struct('name', {'dPp'}, 'idx', {{t,j,k}}); om.add_quad_cost('Cpp', {t,j,k}, [], c, 0, vs); - c = mdi.Delta_T * w * baseMVA * mdi.offer(t).NegativeActiveDeltaPrice(:); + c = w * baseMVA * mdi.offer(t).NegativeActiveDeltaPrice(:); vs = struct('name', {'dPm'}, 'idx', {{t,j,k}}); om.add_quad_cost('Cpm', {t,j,k}, [], c, 0, vs); % weighted fixed reserves cost if mdi.IncludeFixedReserves - c = mdi.Delta_T * w * mdi.FixedReserves(t,j,k).cost(r.igr) * baseMVA; + c = w * mdi.FixedReserves(t,j,k).cost(r.igr) * baseMVA; vs = struct('name', {'R'}, 'idx', {{t,j,k}}); om.add_quad_cost('Rcost', {t,j,k}, [], c, 0, vs); end diff --git a/lib/mostver.m b/lib/mostver.m index 6c2257c..2ea3fd3 100644 --- a/lib/mostver.m +++ b/lib/mostver.m @@ -19,7 +19,7 @@ v = struct( 'Name', 'MOST', ... 'Version', '1.2+', ... 'Release', '', ... - 'Date', '24-Oct-2023' ); + 'Date', '25-Oct-2023' ); if nargout > 0 if nargin > 0 rv = v;