diff --git a/CHANGES.md b/CHANGES.md index 8862de7..16f8d41 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -5,13 +5,20 @@ Change history for MOST since 1.2 --------- +#### 10/24/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. + *Thanks to Stefano Nicolin.* + #### 10/4/23 - Fix [issue #37][9] which caused a fatal error in storage input checks with multiple storage units under some circumstances. *Thanks to Keir Steegstra.* #### 2/3/23 - - Remove extra column in ExpectedRampCost and ignore for single period. + - Remove extra column in mdo.results.ExpectedRampCost and ignore for + single period. Version 1.2 - *Dec 13, 2022* @@ -310,3 +317,4 @@ Version 1.0 - *Jun 1, 2016* [7]: https://arxiv.org/abs/2204.08140 [8]: https://github.com/MATPOWER/most/issues/29 [9]: https://github.com/MATPOWER/most/issues/37 +[10]: https://github.com/MATPOWER/most/issues/39 diff --git a/docs/src/MOST-manual/MOST-manual.tex b/docs/src/MOST-manual/MOST-manual.tex index 7425296..28d6ff9 100644 --- a/docs/src/MOST-manual/MOST-manual.tex +++ b/docs/src/MOST-manual/MOST-manual.tex @@ -858,30 +858,30 @@ \subsubsection{Objective Function} \item[--] expected cost of active power dispatch and redispatch \end{itemize} \begin{align} -f_p(p,p_+,p_-) &= \sum_{t\in T} \sum_{j\in J^t} \sum_{k\in K^{tj}} \psi_\alpha^{tjk} \sum_{i\in I^{tjk}} +f_p(p,p_+,p_-) &= \Delta \cdot \sum_{t\in T} \sum_{j\in J^t} \sum_{k\in K^{tj}} \psi_\alpha^{tjk} \sum_{i\in I^{tjk}} \Bigl[\widetilde{C}_P^{ti}(p^{tijk}) + C_{P+}^{ti}(p_+^{tijk}) + C_{P-}^{ti}(p_-^{tijk}) \Bigr] \label{eq:most_energy_cost} \end{align} \begin{itemize} \item[--] cost of zonal reserves\footnotemark[\value{footnote}] \begin{equation} -f_z(r_z) = \sum_{t\in T} \sum_{j\in J^t} \sum_{k\in K^{tj}} \psi_\alpha^{tjk} \sum_{i\in I^{tjk}} C_z^{ti}(r_z^{tijk}) +f_z(r_z) = \Delta \cdot \sum_{t\in T} \sum_{j\in J^t} \sum_{k\in K^{tj}} \psi_\alpha^{tjk} \sum_{i\in I^{tjk}} C_z^{ti}(r_z^{tijk}) \label{eq:most_zres_cost} \end{equation} \item[--] cost of endogenous contingency reserves\footnotemark[\value{footnote}] \begin{equation} -f_r(r_+, r_-) = \sum_{t\in T} \gamma^t \sum_{i\in I^t} \left[C_{R+}^{ti}(r_+^{ti}) + C_{R-}^{ti}(r_-^{ti}) \right] +f_r(r_+, r_-) = \Delta \cdot \sum_{t\in T} \gamma^t \sum_{i\in I^t} \left[C_{R+}^{ti}(r_+^{ti}) + C_{R-}^{ti}(r_-^{ti}) \right] \label{eq:most_cres_cost} \end{equation} \item[--] expected cost of load-following ramping (wear and tear) \begin{equation} -f_\delta(p) = \sum_{t\in T} \! \gamma^t \!\!\!\!\! \sum_{\begin{aligned}\scriptstyle j_1 &\scriptstyle \in J^{t-1}\\[-6pt] \scriptstyle j_2 &\scriptstyle \in J^t\end{aligned}} \!\!\!\!\! +f_\delta(p) = \Delta \cdot \sum_{t\in T} \! \gamma^t \!\!\!\!\! \sum_{\begin{aligned}\scriptstyle j_1 &\scriptstyle \in J^{t-1}\\[-6pt] \scriptstyle j_2 &\scriptstyle \in J^t\end{aligned}} \!\!\!\!\! \phi^{t j_2 j_1} \!\!\!\! \sum_{i\in I^{tj_2 0}} \!\!\! C_\delta^i(p^{tij_20} - p^{(t-1)ij_10}) \label{eq:rampcost} \end{equation} \item[--] cost of load-following ramp reserves \begin{equation} -f_{\rm lf}(\delta_+, \delta_-) = \sum_{t\in T} \gamma^t \sum_{i\in I^t} \left[C_{\delta+}^{ti}(\delta_+^{ti}) + C_{\delta-}^{ti}(\delta_-^{ti}) \right] +f_{\rm lf}(\delta_+, \delta_-) = \Delta \cdot \sum_{t\in T} \gamma^t \sum_{i\in I^t} \left[C_{\delta+}^{ti}(\delta_+^{ti}) + C_{\delta-}^{ti}(\delta_-^{ti}) \right] \label{eq:most_rampres_cost} \end{equation} \item[--] cost of initial stored energy and value (since it is negative) of expected leftover stored energy in terminal states @@ -890,7 +890,7 @@ \subsubsection{Objective Function} \end{equation} \item[--] no load, startup and shutdown costs \begin{equation} -f_{\rm uc}(u,v,w) = \sum_{t\in T} \gamma^t \sum_{i\in I^t} \Bigl( C_P^{ti}(0) u^{ti} + C_v^{ti} v^{ti} + C_w^{ti} w^{ti} \Bigr) +f_{\rm uc}(u,v,w) = \sum_{t\in T} \gamma^t \sum_{i\in I^t} \Bigl( \Delta \cdot C_P^{ti}(0) u^{ti} + C_v^{ti} v^{ti} + C_w^{ti} w^{ti} \Bigr) \end{equation} \end{itemize} @@ -3441,6 +3441,7 @@ \subsubsection*{Changes} \subsubsection*{Bugs Fixed} \begin{itemize} +\item Fix issue \#39 in which the value of \code{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. \emph{Thanks to Stefano Nicolin.} \item Fix issue \#37 which caused a fatal error in storage input checks with multiple storage units under some circumstances. \emph{Thanks to Keir Steegstra.} \end{itemize} diff --git a/lib/most.m b/lib/most.m index 0cba263..edc77ec 100644 --- a/lib/most.m +++ b/lib/most.m @@ -56,9 +56,8 @@ % MDO MOST data structure, output % (see MOST User's Manual for details) - % MOST -% Copyright (c) 2010-2022, Power Systems Engineering Research Center (PSERC) +% Copyright (c) 2010-2023, Power Systems Engineering Research Center (PSERC) % by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia % and Ray Zimmerman, PSERC Cornell % @@ -82,7 +81,7 @@ fprintf( ' ----- Built on MATPOWER -----\n'); fprintf( ' by Carlos E. Murillo-Sanchez, Universidad Nacional de Colombia--Manizales\n'); fprintf( ' and Ray D. Zimmerman, Cornell University\n'); - fprintf( ' (c) 2012-2022 Power Systems Engineering Research Center (PSERC) \n'); + fprintf( ' (c) 2012-2023 Power Systems Engineering Research Center (PSERC) \n'); fprintf( '=============================================================================\n'); end @@ -509,7 +508,7 @@ for k = 1:mdi.idx.nc(t,j)+1 mpc = mdi.flow(t,j,k).mpc; c00tjk = totcost(mpc.gencost, zeros(ng,1)); - mdi.UC.c00(:, t) = mdi.UC.c00(:, t) + mdi.CostWeightsAdj(k, j, t) * c00tjk; + mdi.UC.c00(:, t) = mdi.UC.c00(:, t) + mdi.Delta_T * mdi.CostWeightsAdj(k, j, t) * c00tjk; c0col = COST + mpc.gencost(:,NCOST) - 1; ipoly = find(mpc.gencost(:, MODEL) == POLYNOMIAL); ipwl = find(mpc.gencost(:, MODEL) == PW_LINEAR); @@ -1724,10 +1723,10 @@ 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 = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,1), 0, ng, ng); - c = -w * baseMVA * mdi.RampWearCostCoeff(:,1) .* mdi.InitialPg; + 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; vs = struct('name', {'Pg'}, 'idx', {{1,j,1}}); - k0 = w * 0.5 * mdi.RampWearCostCoeff(:,1)' * mdi.InitialPg.^2; om.add_quad_cost('RampWear', {1,j,1}, Q, c, k0, vs); end % Then the remaining periods @@ -1735,7 +1734,7 @@ 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 = w * baseMVA^2 * mdi.RampWearCostCoeff(:,t); + h = mdi.Delta_T * w * baseMVA^2 * mdi.RampWearCostCoeff(:,t); 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); @@ -1751,10 +1750,10 @@ if ~mdi.OpenEnded for j = 1:mdi.idx.nj(nt) w = mdi.tstep(nt+1).TransMat(1, j) * mdi.CostWeights(1, j, nt); - Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,nt+1), 0, ng, ng); - c = -w * baseMVA * mdi.RampWearCostCoeff(:,nt+1) .* mdi.TerminalPg; + 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; vs = struct('name', {'Pg'}, 'idx', {{nt,j,1}}); - k0 = w * 0.5 * mdi.RampWearCostCoeff(:,nt+1)' * mdi.TerminalPg.^2; om.add_quad_cost('RampWear', {nt+1,j,1}, Q, c, k0, vs); end end @@ -1785,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, 2 * w * baseMVA^2*gc(ipol, COST), ng, ng); + Q = sparse(ipol, ipol, mdi.Delta_T * 2 * w * baseMVA^2*gc(ipol, COST), ng, ng); else Q = sparse(ng,ng); end c = zeros(ng, 1); if ncost >= 2 - c(ipol) = w * baseMVA*gc(ipol, COST+ncost-2); + c(ipol) = mdi.Delta_T * w * baseMVA*gc(ipol, COST+ncost-2); end - k0 = w * sum(gc(ipol, COST+ncost-1)); + k0 = mdi.Delta_T * w * sum(gc(ipol, COST+ncost-1)); else %% non-uniform order of polynomials %% use a loop Q = sparse(ng,ng); @@ -1803,12 +1802,12 @@ if ncost > 3 error('most: polynomial generator costs of order higher than quadratic not supported'); elseif ncost == 3 - Q(i,i) = 2 * w * baseMVA^2*gc(i, COST); + Q(i,i) = mdi.Delta_T * 2 * w * baseMVA^2*gc(i, COST); end if ncost >= 2 - c(i) = w * baseMVA*gc(i, COST+ncost-2); + c(i) = mdi.Delta_T * w * baseMVA*gc(i, COST+ncost-2); end - k0 = w * gc(i, COST+ncost-1); + k0 = mdi.Delta_T * w * gc(i, COST+ncost-1); end end vs = struct('name', {'Pg'}, 'idx', {{t,j,k}}); @@ -1818,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 = w * ones(mdi.idx.ny(t,j,k),1); + c = mdi.Delta_T * 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 = w * baseMVA * mdi.offer(t).PositiveActiveDeltaPrice(:); + c = mdi.Delta_T * 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 = w * baseMVA * mdi.offer(t).NegativeActiveDeltaPrice(:); + c = mdi.Delta_T * 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 = w * mdi.FixedReserves(t,j,k).cost(r.igr) * baseMVA; + c = mdi.Delta_T * 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 @@ -1841,10 +1840,10 @@ end % contingency reserve costs - c = baseMVA * mdi.StepProb(t) * mdi.offer(t).PositiveActiveReservePrice(:); + c = mdi.Delta_T * baseMVA * mdi.StepProb(t) * mdi.offer(t).PositiveActiveReservePrice(:); vs = struct('name', {'Rpp'}, 'idx', {{t}}); om.add_quad_cost('Crpp', {t}, [], c, 0, vs); - c = baseMVA * mdi.StepProb(t) * mdi.offer(t).NegativeActiveReservePrice(:); + c = mdi.Delta_T * baseMVA * mdi.StepProb(t) * mdi.offer(t).NegativeActiveReservePrice(:); vs = struct('name', {'Rpm'}, 'idx', {{t}}); om.add_quad_cost('Crpm', {t}, [], c, 0, vs); end @@ -1852,10 +1851,10 @@ om.init_indexed_name('qdc', 'Crrp', {mdi.idx.ntramp}); om.init_indexed_name('qdc', 'Crrm', {mdi.idx.ntramp}); for t = 1:mdi.idx.ntramp - c = baseMVA * mdi.StepProb(t) * mdi.offer(t).PositiveLoadFollowReservePrice(:); + c = mdi.Delta_T * baseMVA * mdi.StepProb(t) * mdi.offer(t).PositiveLoadFollowReservePrice(:); vs = struct('name', {'Rrp'}, 'idx', {{t}}); om.add_quad_cost('Crrp', {t}, [], c, 0, vs); - c = baseMVA * mdi.StepProb(t) * mdi.offer(t).NegativeLoadFollowReservePrice(:); + c = mdi.Delta_T * baseMVA * mdi.StepProb(t) * mdi.offer(t).NegativeLoadFollowReservePrice(:); vs = struct('name', {'Rrm'}, 'idx', {{t}}); om.add_quad_cost('Crrm', {t}, [], c, 0, vs); end diff --git a/lib/mostver.m b/lib/mostver.m index 28b6e67..6c2257c 100644 --- a/lib/mostver.m +++ b/lib/mostver.m @@ -19,7 +19,7 @@ v = struct( 'Name', 'MOST', ... 'Version', '1.2+', ... 'Release', '', ... - 'Date', '30-May-2023' ); + 'Date', '24-Oct-2023' ); if nargout > 0 if nargin > 0 rv = v; diff --git a/lib/t/t_most_uc.m b/lib/t/t_most_uc.m index f077acd..0a8507a 100644 --- a/lib/t/t_most_uc.m +++ b/lib/t/t_most_uc.m @@ -43,7 +43,7 @@ function t_most_uc(quiet, create_plots, create_pdfs, savedir) % fcn = {'gurobi'}; % solvers = {'MOSEK'}; % fcn = {'mosek'}; -ntests = 68; +ntests = 69; t_begin(ntests*length(solvers), quiet); if quiet @@ -235,6 +235,10 @@ function t_most_uc(quiet, create_plots, create_pdfs, savedir) plot_case('+ DC Network', mdo, ms, 500, 100, savedir, pp, fname); end % keyboard; + mdi.Delta_T = 2; + mdo = most(mdi, mpopt); + ms = most_summary(mdo); + t_is(ms.f, 2 * ex.f, 8, [t '(Delta_T = 2) : f']); t = sprintf('%s : + startup/shutdown costs : ', solvers{s}); if mpopt.out.all