Fix #39. Properly account for Delta_T in objective function terms.

Ref: #39

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

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}
 

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,18 +1723,18 @@
   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
   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 = 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,44 +1817,44 @@
         % 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
       end
     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
   % Assign load following ramp reserve costs.  Do first nt periods first
   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

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;

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