Merge branch 'EfficiencyImprovements' into 'master'

Updates related to the loop over group reserve horizon when recording results....

Closes #148

See merge request !18

inc/1b_sets.gms

@@ -219,3 +219,13 @@ alias(group, group_);
 
 
 
+sets 
+    tt_agg_circular(t, t, t) "Alternative aggregation ordering with embedded circulars"
+    startp(t) "Temporary subset for time steps"
+    s_realized(s) "All s among realized sft (redundant if always equivalent to s)"
+    sft_resdgn(restype, up_down, grid, node, s,f,t) "Temporary tuplet for reserves by restypeDirectionGridNode"
+;
+
+
+
+

inc/3b_periodicLoop.gms

@@ -545,6 +545,12 @@ option clear = tt_aggregate;
 tt_aggregate(t_current(t+dt_active(t)), tt(t))
     = yes;
 
+// Make alternative aggregation ordering
+option clear=tt_agg_circular; tt_agg_circular(t, t_+dt_circular(t_), t_) $= tt_aggregate(t, t_);
+$macro tt_aggcircular(t, t_)  tt_agg_circular(t, t_, t__)
+*$macro tt_aggcircular(t, t_) (tt_aggregate(t, t__), t_(t__+dt_circular(t__)))
+
+
 * =============================================================================
 * --- Defining unit aggregations and ramps ------------------------------------
 * =============================================================================

inc/3c_inputsLoop.gms

@@ -329,32 +329,36 @@ loop(cc(counter),
 
     // Select and average time series data matching the intervals
     ts_unit_(unit_timeseries(unit), param_unit, ft(f, tt_interval(t)))
-        = sum(tt_aggregate(t, t_),
-            ts_unit(unit, param_unit, f, t_+dt_circular(t_))
+        = sum(tt_aggcircular(t, t_),
+            ts_unit(unit, param_unit, f, t_)
             )
             / mInterval(mSolve, 'stepsPerInterval', counter);
 $ontext
 * Should these be handled here at all? See above comment
     ts_effUnit_(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, ft(f, tt_interval(t)))
-        = sum(tt_aggregate(t, t_),
-            ts_effUnit(effSelector, unit, effSelector, param_eff, f, t_+dt_circular(t_))
+        = sum(tt_aggcircular(t, t_),
+            ts_effUnit(effSelector, unit, effSelector, param_eff, f, t_)
             )
             / mInterval(mSolve, 'stepsPerInterval', counter);
     ts_effGroupUnit_(effSelector, unit_timeseries(unit), param_eff, ft(f, tt_interval(t)))
-        = sum(tt_aggregate(t, t_),
-            ts_effGroupUnit(effSelector, unit, param_eff, f, t_+dt_circular(t_))
+        = sum(tt_aggcircular(t, t_),
+            ts_effGroupUnit(effSelector, unit, param_eff, f, t_)
             )
             / mInterval(mSolve, 'stepsPerInterval', counter);
 $offtext
-    ts_influx_(gn(grid, node), sft(s, f, tt_interval(t)))
+    ts_influx_(gn, sft(s, f, tt_interval(t)))$gn_scenarios(gn, 'ts_influx')
         = sum(tt_aggregate(t, t_),
-            ts_influx(grid, node,
-                f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_influx'))
-                     + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_influx')),
-                t_+ (+ dt_scenarioOffset(grid, node, 'ts_influx', s)
-                     + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_influx'))))
-            )
-            / mInterval(mSolve, 'stepsPerInterval', counter);
+            ts_influx(gn,
+                f + (  df_realization(f, t)$(not gn_forecasts(gn, 'ts_influx'))
+                     + df_scenario(f, t)),
+                t_+ (+ dt_scenarioOffset(gn, 'ts_influx', s)))
+            ) / mInterval(mSolve, 'stepsPerInterval', counter);
+    ts_influx_(gn, sft(s, f, tt_interval(t)))$(not gn_scenarios(gn, 'ts_influx'))
+        = sum(tt_aggcircular(t, t_),
+            ts_influx(gn,
+                f + (  df_realization(f, t)$(not gn_forecasts(gn, 'ts_influx'))),
+                t_  )
+            ) / mInterval(mSolve, 'stepsPerInterval', counter);
     ts_cf_(flowNode(flow, node), sft(s, f, tt_interval(t)))
         = sum(tt_aggregate(t, t_),
             ts_cf(flow, node,
@@ -367,10 +371,10 @@ $offtext
     // Reserves relevant only until reserve_length
     ts_reserveDemand_(restypeDirectionGroup(restype, up_down, group), ft(f, tt_interval(t)))
       ${ord(t) <= tSolveFirst + p_groupReserves(group, restype, 'reserve_length')  }
-        = sum(tt_aggregate(t, t_),
+        = sum(tt_aggcircular(t, t_),
             ts_reserveDemand(restype, up_down, group,
                 f + (  df_realization(f, t)${not sum(gnGroup(grid, node, group), gn_forecasts(restype, node, 'ts_reserveDemand'))}
-                     + df_scenario(f, t)${sum(gnGroup(grid, node, group), gn_scenarios(restype, node, 'ts_reserveDemand'))} ),                t_+ dt_circular(t_))
+                     + df_scenario(f, t)${sum(gnGroup(grid, node, group), gn_scenarios(restype, node, 'ts_reserveDemand'))} ), t_)
             )
             / mInterval(mSolve, 'stepsPerInterval', counter);
     ts_node_(gn_state(grid, node), param_gnBoundaryTypes, sft(s, f, tt_interval(t)))
@@ -404,6 +408,7 @@ $offtext
                             + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node'))))
                 )
                 $(sameas(param_gnBoundaryTypes, 'upwardLimit') or upwardSlack(param_gnBoundaryTypes));
+                
     ts_gnn_(gn2n_timeseries(grid, node, node_, param_gnn), ft(f, tt_interval(t)))
         = sum(tt_aggregate(t, t_),
             ts_gnn(grid, node, node_, param_gnn, f, t_+dt_circular(t_))
@@ -414,17 +419,17 @@ $offtext
     ts_vomCost_(gnu(grid, node, unit), tt_interval(t))
         = + p_gnu(grid, node, unit, 'vomCosts')
           // input node cost
-          + p_price(node, 'price')${un_commodity_in(unit, node) and p_price(node, 'useConstant')}
-          + sum(tt_aggregate(t, t_)${un_commodity_in(unit, node) and p_price(node, 'useTimeSeries')},
-              + ts_price(node, t_+dt_circular(t_))
-              )
-              / mInterval(mSolve, 'stepsPerInterval', counter)
+          + (
+             + p_price(node, 'price')$p_price(node, 'useConstant')
+             + sum(tt_aggcircular(t, t_), ts_price(node, t_))$p_price(node, 'useTimeSeries')
+                 / mInterval(mSolve, 'stepsPerInterval', counter)
+            )$un_commodity_in(unit, node) 
           // output node cost (if price > 0 --> ts_vomCost_ < 0, i.e. considered as revenue)
-          - p_price(node, 'price')${un_commodity_out(unit, node) and p_price(node, 'useConstant')}
-          - sum(tt_aggregate(t, t_)${un_commodity_out(unit, node) and p_price(node, 'useTimeSeries')},
-              + ts_price(node, t_+dt_circular(t_))
-              )
-              / mInterval(mSolve, 'stepsPerInterval', counter)
+          - (
+             + p_price(node, 'price')$p_price(node, 'useConstant')
+             + sum(tt_aggcircular(t, t_), ts_price(node, t_))$p_price(node, 'useTimeSeries')
+                 / mInterval(mSolve, 'stepsPerInterval', counter)
+            )$un_commodity_out(unit, node)
           // emission cost
           + sum(emission$p_unitEmissionCost(unit, node, emission), // Emission taxes
               + p_unitEmissionCost(unit, node, emission)
@@ -444,9 +449,9 @@ $offtext
             + p_unStartup(unit, node, starttype) // MWh/start-up
               * [
                   + p_price(node, 'price')$p_price(node, 'useConstant') // CUR/MWh
-                  + sum(tt_aggregate(t, t_)$p_price(node, 'useTimeseries'),
-                      + ts_price(node, t_+dt_circular(t_)) // CUR/MWh
-                    )
+                  + sum(tt_aggcircular(t, t_),
+                      + ts_price(node, t_) // CUR/MWh
+                    )$p_price(node, 'useTimeseries')
                     / mInterval(mSolve, 'stepsPerInterval', counter)
                 ] // END * p_uStartup
           ) // END sum(node)

inc/4a_outputVariant.gms

@@ -31,110 +31,112 @@ if(tSolveFirst >= mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializa
               + v_online_MIP.l(unit, s, f_solve, t)$unit_online_MIP(unit);
     );
 );
+
+// Improve performance & readibility by using a few helper sets
+option clear=startp, clear=sft_resdgn, s_realized < sft_realized;
+startp(t)$(ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod'))=yes;
+
 // Realized state history
-loop(sft_realized(s, f, t),
-    r_state(gn_state(grid, node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
+loop(ms(mSolve, s_realized(s)),
+    r_state(gn_state(grid, node), f, startp(t))$sft_realized(s, f, t)
         = v_state.l(grid, node, s, f, t);
 
     // Realized state history - initial state values in samples
-    r_state(gn_state(grid, node), f_solve(f), t+dt(t))${   ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
-                                                           and sum(ms(mSolve, s), mst_start(mSolve, s, t))
-                                                       }
-        = v_state.l(grid, node, s, f, t+dt(t))
+    r_state(gn_state(grid, node), f_solve(f), t_(t+dt(t)))$(mst_start(ms, t)$sft_realized(s, f, t)$startp(t))
+        = v_state.l(grid, node, s, f, t_)
     ;
     // Realized unit online history
-    r_online(uft_online(unit, f, t))$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
-        = v_online_LP.l(unit, s, f, t)${ uft_onlineLP(unit, f, t)    }
-            + v_online_MIP.l(unit, s, f, t)${  uft_onlineMIP(unit, f, t)   }
+    r_online(uft_online(unit, f, startp(t)))$sft_realized(s, f, t)
+        = v_online_LP.l(unit, s, f, t)$uft_onlineLP(unit, f, t)
+            + v_online_MIP.l(unit, s, f, t)$uft_onlineMIP(unit, f, t)
     ;
     // Unit startup and shutdown history
-    r_startup(unit, starttype, f, t)${ uft_online(unit, f, t) and [ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')] }
-        = v_startup_LP.l(unit, starttype, s, f, t)${ uft_onlineLP(unit, f, t) }
-            + v_startup_MIP.l(unit, starttype, s, f, t)${ uft_onlineMIP(unit, f, t) }
+    r_startup(unit, starttype, f, startp(t))$(uft_online(unit, f, t)$sft_realized(s, f, t))
+        = v_startup_LP.l(unit, starttype, s, f, t)$uft_onlineLP(unit, f, t)
+            + v_startup_MIP.l(unit, starttype, s, f, t)$uft_onlineMIP(unit, f, t)
     ;
-    r_shutdown(uft_online(unit, f, t))$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
-        = v_shutdown_LP.l(unit, s, f, t)${ uft_onlineLP(unit, f, t) }
-            + v_shutdown_MIP.l(unit, s, f, t)${ uft_onlineMIP(unit, f, t) }
+    r_shutdown(uft_online(unit, f, startp(t)))$sft_realized(s, f, t)
+        = v_shutdown_LP.l(unit, s, f, t)$uft_onlineLP(unit, f, t)
+            + v_shutdown_MIP.l(unit, s, f, t)$uft_onlineMIP(unit, f, t)
     ;
 );
 
 * --- Reserve results ---------------------------------------------------------
 
 // Loop over reserve horizon, as the reserve variables use a different ft-structure due to commitment
-loop((restypeDirectionGridNode(restype, up_down, grid, node), sft(s, f, t))
-    ${  ord(t) <= tSolveFirst + p_gnReserves(grid, node, restype, 'reserve_length')
-        and ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
-        },
 
+sft_resdgn(restypeDirectionGridNode(restype, up_down, gn), sft(s, f, startp(t)))
+  ${ord(t) <= tSolveFirst + p_gnReserves(gn, restype, 'reserve_length')} = yes;
+
+loop(s,
     // Reserve provisions of units
-    r_reserve(gnuRescapable(restype, up_down, grid, node, unit), f+df_reserves(grid, node, restype, f, t), t)
-        ${  not [   restypeReleasedForRealization(restype)
-                    and sft_realized(s, f+df_reserves(grid, node, restype, f, t), t)
-                    ]
-            }
-        = + v_reserve.l(restype, up_down, grid, node, unit, s, f+df_reserves(grid, node, restype, f, t), t)
-          + sum(restype_$p_gnuRes2Res(grid, node, unit, restype_, up_down, restype),
-              + v_reserve.l(restype_, up_down, grid, node, unit, s, f+df_reserves(grid, node, restype_, f, t), t)
-                  * p_gnuRes2Res(grid, node, unit, restype_, up_down, restype)
+    r_reserve(gnuRescapable(restype, up_down, gn, unit), f_(f+df_reserves(gn, restype, f, t)), t)
+        ${ (not sft_realized(s, f_, t)$restypeReleasedForRealization(restype))$sft_resdgn(restype,up_down,gn,s,f,t) }
+        = + v_reserve.l(restype, up_down, gn, unit, s, f_, t)
+          + sum(restype_$p_gnuRes2Res(gn, unit, restype_, up_down, restype),
+              + v_reserve.l(restype_, up_down, gn, unit, s, f+df_reserves(gn, restype_, f, t), t)
+                  * p_gnuRes2Res(gn, unit, restype_, up_down, restype)
             );
 
+    // Reserve transfer capacity for links defined out from this node
+    r_resTransferRightward(restype, up_down, gn2n_directional(gn, to_node), f_(f+df_reserves(gn, restype, f, t)), t)
+        ${ restypeDirectionGridNodeNode(restype, up_down, gn, to_node)$sft_resdgn(restype,up_down,gn,s,f,t) }
+        = v_resTransferRightward.l(restype, up_down, gn, to_node, s, f_, t);
+
+    r_resTransferLeftward(restype, up_down, gn2n_directional(gn(grid, node), to_node), f_(f+df_reserves(grid, to_node, restype, f, t)), t)
+        ${ restypeDirectionGridNodeNode(restype, up_down, grid, to_node, node)$sft_resdgn(restype,up_down,gn,s,f,t) }
+        = v_resTransferLeftward.l(restype, up_down, gn, to_node, s, f_, t);
+
+); // END loop(restypeDirectionNode, sft)
+
+// Loop over group reserve horizon
+loop((restypeDirectionGroup(restype, up_down, group), sft(s, f, startp(t)))
+    ${ord(t) <= tSolveFirst + p_groupReserves(group, restype, 'reserve_length')},
+
     // Reserve requirement due to N-1 reserve constraint
-    r_resDemandLargestInfeedUnit(restypeDirectionGroup(restype, 'up', group), f+df_reserves(grid, node, restype, f, t), t)
-        ${ sum((gnGroup(grid, node, group), unit_fail), p_gnuReserves(grid, node, unit_fail, restype, 'portion_of_infeed_to_reserve')) } // Calculate only for groups with units that can fail.
-        = smax(unit_fail(unit_)${ sum(gnGroup(grid, node, group), p_gnuReserves(grid, node, unit_, restype, 'portion_of_infeed_to_reserve')) },
-            + v_gen.l(grid, node, unit_, s, f, t)
-                * p_gnuReserves(grid, node, unit_, restype, 'portion_of_infeed_to_reserve')
+    r_resDemandLargestInfeedUnit(restype, 'up', group, f_(f+df_reservesGroup(group, restype, f, t)), t)
+        ${ sum((gnGroup(gn, group),unit_fail)$p_gnuReserves(gn, unit_fail, restype, 'portion_of_infeed_to_reserve'),1) } // Calculate only for groups with units that can fail.
+        = smax((gnGroup(gn, group),unit_fail)$p_gnuReserves(gn, unit_fail, restype, 'portion_of_infeed_to_reserve'),
+            + v_gen.l(gn, unit_fail, s, f, t)
+                * p_gnuReserves(gn, unit_fail, restype, 'portion_of_infeed_to_reserve')
             ) // END smax(unit_fail)
         ;
-
-    // Reserve transfer capacity for links defined out from this node
-    r_resTransferRightward(restype, up_down, grid, node, to_node, f+df_reserves(grid, node, restype, f, t), t)
-        ${  gn2n_directional(grid, node, to_node)
-            and restypeDirectionGridNodeNode(restype, up_down, grid, node, to_node)
-            }
-        = v_resTransferRightward.l(restype, up_down, grid, node, to_node, s, f+df_reserves(grid, node, restype, f, t), t);
-
-    r_resTransferLeftward(restype, up_down, grid, node, to_node, f+df_reserves(grid, to_node, restype, f, t), t)
-        ${  gn2n_directional(grid, node, to_node)
-            and restypeDirectionGridNodeNode(restype, up_down, grid, to_node, node)
-            }
-        = v_resTransferLeftward.l(restype, up_down, grid, node, to_node, s, f+df_reserves(grid, to_node, restype, f, t), t);
-
     // Dummy reserve demand changes
-    r_qResDemand(restype, up_down, group, f+df_reservesGroup(group, restype, f, t), t)
-        = vq_resDemand.l(restype, up_down, group, s, f+df_reservesGroup(group, restype, f, t), t);
+    r_qResDemand(restype, up_down, group, f_(f+df_reservesGroup(group, restype, f, t)), t)
+        = vq_resDemand.l(restype, up_down, group, s, f_, t);
 
-    r_qResMissing(restype, up_down, group, f+df_reservesGroup(group, restype, f, t), t)
-        = vq_resMissing.l(restype, up_down, group, s, f+df_reservesGroup(group, restype, f, t), t);
+    r_qResMissing(restype, up_down, group, f_(f+df_reservesGroup(group, restype, f, t)), t)
+        = vq_resMissing.l(restype, up_down, group, s, f_, t);
 
-); // END loop(restypeDirectionNode, sft)
+); // END loop(restypeDirectionGroup, sft)
 
 * --- Interesting results -----------------------------------------------------
 
-loop(sft_realized(s, f, t),
+loop(s_realized(s),
+
     // Unit generation and consumption
-    r_gen(gnuft(grid, node, unit, f, t))$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
+    r_gen(gnuft(grid, node, unit, f, startp(t)))$sft_realized(s, f, t)
         = v_gen.l(grid, node, unit, s, f, t)
     ;
     // Fuel use of units
-    r_fuelUse(node, uft(unit_commodity, f, t))$[un_commodity(unit_commodity, node) and ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
-        = - sum(grid$gnu(grid, node, unit_commodity), v_gen.l(grid, node, unit_commodity, s, f, t))
+    r_fuelUse(node, uft(unit_commodity, f, startp(t)))$(un_commodity(unit_commodity, node)$sft_realized(s, f, t))
+        = - sum(gnu(grid, node, unit_commodity), v_gen.l(grid, node, unit_commodity, s, f, t))
     ;
     // Transfer of energy between nodes
-    r_transfer(gn2n(grid, from_node, to_node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
+    r_transfer(gn2n(grid, from_node, to_node), f, startp(t))$sft_realized(s, f, t)
         = v_transfer.l(grid, from_node, to_node, s, f, t)
     ;
     // Transfer of energy from first node to second node
-    r_transferRightward(gn2n_directional(grid, from_node, to_node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
+    r_transferRightward(gn2n_directional(grid, from_node, to_node), f, startp(t))$sft_realized(s, f, t)
         = v_transferRightward.l(grid, from_node, to_node, s, f, t)
     ;
     // Transfer of energy from second node to first node
-    r_transferLeftward(gn2n_directional(grid, to_node, from_node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
+    r_transferLeftward(gn2n_directional(grid, to_node, from_node), f, startp(t))$sft_realized(s, f, t)
         = v_transferLeftward.l(grid, to_node, from_node, s, f, t)
     ;
     // Energy spilled from nodes
-    r_spill(gn(grid, node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
-        = v_spill.l(grid, node, s, f, t)
+    r_spill(gn, f, startp(t))$sft_realized(s, f, t)
+        = v_spill.l(gn, s, f, t)
     ;
 );
 
@@ -144,17 +146,17 @@ r_totalObj(tSolve)
 ;
 
 // q_balance marginal values
-loop(sft_realized(s, f, t),
-    r_balanceMarginal(gn(grid, node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
-        = q_balance.m(grid, node, mSolve, s, f, t)
+loop(s_realized(s),
+    r_balanceMarginal(gn, f, startp(t))$sft_realized(s, f, t)
+        = q_balance.m(gn, mSolve, s, f, t)
     ;
     // q_resDemand marginal values
-    r_resDemandMarginal(restypeDirectionGroup(restype, up_down, group), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
+    r_resDemandMarginal(restypeDirectionGroup(restype, up_down, group), f, startp(t))$sft_realized(s, f, t)
         = q_resDemand.m(restype, up_down, group, s, f, t)
     ;
     // v_stateSlack values for calculation of realized costs later on
-    r_stateSlack(gn_stateSlack(grid, node), slack, f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
-        = v_stateSlack.l(grid, node, slack, s, f, t)
+    r_stateSlack(gn_stateSlack(gn), slack, f, startp(t))$sft_realized(s, f, t)
+        = v_stateSlack.l(gn, slack, s, f, t)
     ;
 );
 // Unit investments
@@ -175,16 +177,12 @@ r_investTransfer(grid, node, node_, t_invest(t))${ p_gnn(grid, node, node_, 'tra
 * --- Feasibility results -----------------------------------------------------
 loop(sft_realized(s, f, t),
 // Dummy generation & consumption
-r_qGen(inc_dec, gn(grid, node), f, t)
-    ${  ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
-        }
-    = vq_gen.l(inc_dec, grid, node, s, f, t)
+r_qGen(inc_dec, gn, f, startp(t))
+    = vq_gen.l(inc_dec, gn, s, f, t)
 ;
 // Dummy capacity
-r_qCapacity(gn(grid, node), f, t)
-    ${  ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
-        }
-    = vq_capacity.l(grid, node, s, f, t)
+r_qCapacity(gn, f, startp(t))
+    = vq_capacity.l(gn, s, f, t)
 ;
 );
 

inc/4b_outputInvariant.gms

@@ -15,6 +15,10 @@ You should have received a copy of the GNU Lesser General Public License
 along with Backbone.  If not, see <http://www.gnu.org/licenses/>.
 $offtext
 
+
+// for performance, get rid of any zeros in r_gen (r_gen was already missing many zero values anyway)
+r_gen(gnu, f, t)$((r_gen(gnu, f, t)=0)$r_gen(gnu, f, t))=0;
+
 * =============================================================================
 * --- Time Step Dependent Results ---------------------------------------------
 * =============================================================================
@@ -22,17 +26,18 @@ $offtext
 // Need to loop over the model dimension, as this file is no longer contained in the modelSolves loop...
 loop(m,
 
+    option clear=startp; startp(t)$(ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod'))=yes;
 * --- Realized Individual Costs ----------------------------------------------
 
     // Variable O&M costs
-    r_gnuVOMCost(gnu(grid, node, unit), ft_realizedNoReset(f,t))$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+    r_gnuVOMCost(gnu(grid, node, unit), ft_realizedNoReset(f,startp(t)))
         = 1e-6 // Scaling to MEUR
             * p_stepLengthNoReset(m, f, t)
             * abs(r_gen(grid, node, unit, f, t))
             * p_gnu(grid, node, unit, 'vomCosts');
 
     // Fuel and emission costs during normal operation
-    r_uFuelEmissionCost(commodity, unit_commodity(unit), ft_realizedNoReset(f,t))${ un_commodity(unit, commodity) and [ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]}
+    r_uFuelEmissionCost(commodity, unit_commodity(unit), ft_realizedNoReset(f,startp(t)))$un_commodity(unit, commodity)
         = 1e-6 // Scaling to MEUR
             * p_stepLengthNoReset(m, f, t)
             * r_fuelUse(commodity, unit, f, t)
@@ -43,7 +48,7 @@ loop(m,
               ];
 
     // Unit startup costs
-    r_uStartupCost(unit, ft_realizedNoReset(f,t))${sum(starttype, unitStarttype(unit, starttype)) and [ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]}
+    r_uStartupCost(unit, ft_realizedNoReset(f,startp(t)))$sum(starttype, unitStarttype(unit, starttype))
         = 1e-6 // Scaling to MEUR
             * sum(unitStarttype(unit, starttype),
                 + r_startup(unit, starttype, f, t)
@@ -65,7 +70,7 @@ loop(m,
               ); // END sum(starttype)
 
     //Variable Trnasfer Costs
-    r_gnnVariableTransCost(gn2n_directional(grid, node_, node), ft_realizedNoReset(f,t))$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+    r_gnnVariableTransCost(gn2n_directional(grid, node_, node), ft_realizedNoReset(f,startp(t)))
         = 1e-6 // Scaling to MEUR
             * p_stepLengthNoReset(m, f, t)
                     *[+ p_gnn(grid, node, node_, 'variableTransCost')
@@ -74,7 +79,7 @@ loop(m,
                     * r_transferRightward(grid, node_, node, f, t)];
 
     // Node state slack costs
-    r_gnStateSlackCost(gn_stateSlack(grid, node), ft_realizedNoReset(f,t))$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+    r_gnStateSlackCost(gn_stateSlack(grid, node), ft_realizedNoReset(f,startp(t)))
         = 1e-6 // Scaling to MEUR
             * p_stepLengthNoReset(m, f, t)
             * sum(slack${ p_gnBoundaryPropertiesForStates(grid, node, slack, 'slackCost') },
@@ -106,35 +111,35 @@ loop(m,
 
     // Total VOM costs
     r_gnuTotalVOMCost(gnu_output(grid, node, unit))
-        = sum(ft_realizedNoReset(f,t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f,startp(t)),
             + r_gnuVOMCost(grid, node, unit, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s) * p_s_discountFactor(s))
             );
 
     // Total Variable Transfer costs
     r_gnnTotalVariableTransCost(gn2n_directional(grid, node_, node))
-        = sum(ft_realizedNoReset(f,t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f,startp(t)),
             + r_gnnVariableTransCost(grid, node_, node, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s) * p_s_discountFactor(s))
             );
 
     // Total fuel & emission costs
     r_uTotalFuelEmissionCost(commodity, unit)$un_commodity(unit, commodity)
-        = sum(ft_realizedNoReset(f,t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f,startp(t)),
             + r_uFuelEmissionCost(commodity, unit, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s) * p_s_discountFactor(s))
             );
 
     // Total unit startup costs
-    r_uTotalStartupCost(unit)${ sum(starttype, unitStarttype(unit, starttype)) }
-        = sum(ft_realizedNoReset(f,t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+    r_uTotalStartupCost(unit)$sum(starttype, unitStarttype(unit, starttype))
+        = sum(ft_realizedNoReset(f,startp(t)),
             + r_uStartupCost(unit, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s) * p_s_discountFactor(s))
             );
 
     // Total state variable slack costs
     r_gnTotalStateSlackCost(gn_stateSlack(grid, node))
-        = sum(ft_realizedNoReset(f,t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f,startp(t)),
             + r_gnStateSlackCost(grid, node, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s) * p_s_discountFactor(s))
             );
@@ -185,19 +190,17 @@ loop(m,
 * --- Realized Nodal System Costs ---------------------------------------------
 
     // Total realized gn operating costs
-    r_gnRealizedOperatingCost(gn(grid, node), ft_realizedNoReset(f, t))$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+    r_gnRealizedOperatingCost(gn(grid, node), ft_realizedNoReset(f, startp(t)))
         = + sum(gnu(grid, node, unit),
               // VOM costs
               + r_gnuVOMCost(grid, node, unit, f, t)
             )
           // Allocate fuel and startup costs on energy basis, but for output nodes only
-          + sum(unit$gnu(grid, node, unit),
-              + [
-                  + abs{r_gen(grid, node, unit, f, t)}  // abs is due to potential negative outputs like energy from a cooling unit. It's the energy contribution that matters, not direction.
-                      / sum(gnu_output(grid_output, node_output, unit),
-                          + abs{r_gen(grid_output, node_output, unit, f, t)}
-                        ) // END sum(gnu_output)
-                ]$(gnu_output(grid, node, unit) and abs{r_gen(grid, node, unit, f, t)} > eps)
+          + sum(unit$(r_gen(grid, node, unit, f, t)$gnu_output(grid, node, unit)),
+              + abs{r_gen(grid, node, unit, f, t)}  // abs is due to potential negative outputs like energy from a cooling unit. It's the energy contribution that matters, not direction.
+                   / sum(gnu_output(grid_output, node_output, unit),
+                       + abs{r_gen(grid_output, node_output, unit, f, t)}
+                     ) // END sum(gnu_output)
                 *
                 {
                   + sum(un_commodity(unit, commodity), r_uFuelEmissionCost(commodity, unit, f, t))
@@ -214,7 +217,7 @@ loop(m,
 * --- Realized Nodal Energy Consumption ---------------------------------------
 // !!! NOTE !!! This is a bit of an approximation at the moment !!!!!!!!!!!!!!!
 
-    r_gnConsumption(gn(grid, node), ft_realizedNoReset(f, t))$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+    r_gnConsumption(gn(grid, node), ft_realizedNoReset(f, startp(t)))
         = p_stepLengthNoReset(m, f, t)
             * [
                 + min(ts_influx(grid, node, f, t), 0) // Not necessarily a good idea, as ts_influx contains energy gains as well...
@@ -227,7 +230,7 @@ loop(m,
 
     // Total energy generation
     r_gnuTotalGen(gnu_output(grid, node, unit))
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_gen(grid, node, unit, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -235,7 +238,7 @@ loop(m,
 
     r_gen_gnUnittype(gn(grid, node), unittype)$sum(unit$unitUnittype(unit, unittype), 1)
       = sum(gnu(grid,node,unit)$unitUnittype(unit, unittype),
-            sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+            sum(ft_realizedNoReset(f, startp(t)),
               + r_gen(grid, node, unit, f, t)
                   * p_stepLengthNoReset(m, f, t)
                   * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -243,8 +246,7 @@ loop(m,
         );
 
     // Energy generation by fuels
-    r_genFuel(gn(grid, node), commodity, ft_realizedNoReset(f, t))$[    sum(gnu_input(grid_, node_, unit)$gnu_output(grid, node, unit), r_gen(grid_, node_, unit, f, t))
-                                                                    and ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+    r_genFuel(gn(grid, node), commodity, ft_realizedNoReset(f, startp(t)))$sum(gnu_input(grid_, commodity, unit)$gnu_output(grid, node, unit),r_gen(grid_, commodity, unit, f, t))
         = sum(gnu_output(grid, node, unit)$sum(gnu_input(grid_, commodity, unit), 1),
             + r_gen(grid, node, unit, f, t)
           );
@@ -258,17 +260,14 @@ loop(m,
             + r_gen(grid, node, unit, f, t));
 
     // Energy generation by fuels
-    r_genUnittype(gn(grid, node), unittype, ft_realizedNoReset(f,t))
-        ${  sum(unit,gnu_output(grid, node, unit))
-            and [ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
-            }
-        = sum(unit${unitUnittype(unit, unittype) and gnu_output(grid, node, unit)},
+    r_genUnittype(gn(grid, node), unittype, ft_realizedNoReset(f,startp(t)))
+        = sum(gnu_output(grid, node, unit)$unitUnittype(unit, unittype),
             + r_gen(grid, node, unit, f, t)
             ); // END sum(unit)
 
     // Total generation on each node by fuels
     r_gnTotalGenFuel(gn(grid, node), commodity)
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_genFuel(grid, node, commodity, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -276,7 +275,7 @@ loop(m,
 
     // Total dummy generation/consumption
     r_gnTotalqGen(inc_dec, gn(grid, node))
-        = sum(ft_realizedNoReset(f,t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f,startp(t)),
             + r_qGen(inc_dec, grid, node, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -285,7 +284,7 @@ loop(m,
 * --- Emission Results --------------------------------------------------------
 
     // Emissions of units (only for commodities, not including startup fuels)
-    r_emissions(commodity, emission, unit, ft_realizedNoReset(f,t))$ {un_commodity(unit, commodity) and [ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]}
+    r_emissions(commodity, emission, unit, ft_realizedNoReset(f,startp(t)))$un_commodity(unit, commodity)
         =   + p_stepLengthNoReset(m, f, t)
             * r_fuelUse(commodity, unit, f, t)
             * p_nEmission(commodity, emission)
@@ -295,7 +294,7 @@ loop(m,
 
     // Emission sums
     r_nuTotalEmissions (commodity, unit, emission)
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_emissions(commodity, emission, unit, f, t)
                  * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
             )
@@ -317,7 +316,7 @@ loop(m,
 
     // Total sub-unit-hours for units over the simulation
     r_uTotalOnline(unit)
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_online(unit, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -340,7 +339,7 @@ loop(m,
 
     // Total reserve provisions over the simulation
     r_gnuTotalReserve(gnuRescapable(restype, up_down, grid, node, unit))
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_reserve(restype, up_down, grid, node, unit, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -348,7 +347,7 @@ loop(m,
 
     // Total dummy reserve provisions over the simulation
     r_groupTotalqResDemand(restypeDirectionGroup(restype, up_down, group))
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_qResDemand(restype, up_down, group, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -358,7 +357,7 @@ loop(m,
 
     // Total transfer of energy between nodes
     r_gnnTotalTransfer(gn2n(grid, from_node, to_node))
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_transfer(grid, from_node, to_node, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -366,7 +365,7 @@ loop(m,
 
     // Total energy spill from nodes
     r_gnTotalSpill(grid, node_spill(node))
-        = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+        = sum(ft_realizedNoReset(f, startp(t)),
             + r_spill(grid, node, f, t)
                 * p_stepLengthNoReset(m, f, t)
                 * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
@@ -379,11 +378,11 @@ loop(m,
 * --- Scaling Marginal Values to EUR/MWh from MEUR/MWh ------------------------
 
 // Energy balance
-r_balanceMarginal(gn(grid, node), ft_realizedNoReset(f, t))$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+r_balanceMarginal(gn(grid, node), ft_realizedNoReset(f, startp(t)))
     = 1e6 * r_balanceMarginal(grid, node, f, t);
 
 // Reserve balance
-r_resDemandMarginal(restypeDirectionGroup(restype, up_down, group), ft_realizedNoReset(f, t))$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')]
+r_resDemandMarginal(restypeDirectionGroup(restype, up_down, group), ft_realizedNoReset(f, startp(t)))
     = 1e6 * r_resDemandMarginal(restype, up_down, group, f, t);
 
 * --- Total Generation Results ------------------------------------------------
@@ -416,7 +415,7 @@ r_gTotalqGen(inc_dec, grid)
 
 // Total consumption on each gn over the simulation
 r_gnTotalConsumption(gn(grid, node))
-    = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+    = sum(ft_realizedNoReset(f, startp(t)),
         + r_gnConsumption(grid, node, f ,t)
             * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
         );
@@ -460,7 +459,7 @@ r_gnTotalSpillShare(gn(grid, node_spill))${ r_gTotalSpill(grid) > 0 }
 
 // Total realized operating costs on each gn over the simulation
 r_gnTotalRealizedOperatingCost(gn(grid, node))
-    = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+    = sum(ft_realizedNoReset(f, startp(t)),
         + r_gnRealizedOperatingCost(grid, node, f ,t)
             * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s) * p_s_discountFactor(s))
         );
@@ -581,14 +580,14 @@ r_uTotalOnlinePerUnit(unit)${ p_unit(unit, 'unitCount') > 0 }
 
 // Total sub-unit startups over the simulation
 r_uTotalStartup(unit, starttype)
-    = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+    = sum(ft_realizedNoReset(f, startp(t)),
         + r_startup(unit, starttype, f, t)
             * sum(msft_realizedNoReset(m, s, f, t), p_msProbability(m, s) * p_msWeight(m, s))
         ); // END sum(ft_realizedNoReset)
 
 // Total sub-unit shutdowns over the simulation
 r_uTotalShutdown(unit)
-    = sum(ft_realizedNoReset(f, t)$[ord(t) > mSettings(m, 't_start') + mSettings(m, 't_initializationPeriod')],
+    = sum(ft_realizedNoReset(f, startp(t)),