Add new set that defines which nodes and timeseries use forecasts

Set `gn_forecasts(*, node, timeseries)` tells for which grid/flow/restype, node and timeseries combinations use forecasts and by default includes everything.

See: #99

CHANGELOG.md

@@ -7,6 +7,7 @@ All notable changes to this project will be documented in this file.
 - All input files, including *inputData.gdx*, are optional
 - Enabling different combinations of LP and MIP online and invest variables
 - Separate availability parameter for output units in the capacity margin constraint
+- Parameter `gn_forecasts(*, node, timeseries)` to tell which nodes and timeseries use forecasts
 
 ### Changed 
 - Updated tool defintions for Sceleton Titan and Spine Toolbox

defModels/investInit_temp.gms

@@ -82,6 +82,10 @@ if (mType('invest'),
     // Define the number of forecasts used by the model
     mSettings('invest', 'forecasts') = 0;
 
+    // Define which nodes and timeseries use forecasts
+    //Option clear = gn_forecasts;  // By default includes everything, so clear first
+    //gn_forecasts('wind', 'XXX', 'ts_cf') = yes;
+
     // Define forecast properties and features
     mSettings('invest', 't_forecastStart') = 0;                // At which time step the first forecast is available ( 1 = t000001 )
     mSettings('invest', 't_forecastLengthUnchanging') = 0;     // Length of forecasts in time steps - this does not decrease when the solve moves forward (requires forecast data that is longer than the horizon at first)

defModels/scheduleInit_temp.gms

@@ -106,6 +106,10 @@ if (mType('schedule'),
     // Define the number of forecasts used by the model
     mSettings('schedule', 'forecasts') = 3;
 
+    // Define which nodes and timeseries use forecasts
+    //Option clear = gn_forecasts;  // By default includes everything, so clear first
+    //gn_forecasts('wind', 'XXX', 'ts_cf') = yes;
+
     // Define forecast properties and features
     mSettings('schedule', 't_forecastStart') = 1;                  // At which time step the first forecast is available ( 1 = t000001 )
     mSettings('schedule', 't_forecastLengthUnchanging') = 36;      // Length of forecasts in time steps - this does not decrease when the solve moves forward (requires forecast data that is longer than the horizon at first)

inc/1b_sets.gms

@@ -49,7 +49,7 @@ Sets
     unit_investLP(unit) "Units with continuous investments allowed"
     unit_investMIP(unit) "Units with integer investments allowed"
     unit_timeseries(unit) "Units with time series enabled"
-    unit_incHRAdditionalConstraints(unit) "Units that use the two additional incremental heat rate constraints" 
+    unit_incHRAdditionalConstraints(unit) "Units that use the two additional incremental heat rate constraints"
 
 * --- Nodes -------------------------------------------------------------------
     node_spill(node) "Nodes that can spill; used to remove v_spill variables where not relevant"
@@ -128,7 +128,8 @@ Sets
     s_prev(s) "Temporary set for previous sample"
 $if defined scenario
     s_scenario(s, scenario) "Which samples belong to which scenarios"
-    gn_scenarios(*, node, *) "Which grid/flow, node and timeseries/param have data for long-term scenarios"
+    gn_forecasts(*, node, timeseries) "Which grid/flow, node and timeseries use short-term forecasts"
+    gn_scenarios(*, node, timeseries) "Which grid/flow, node and timeseries have data for long-term scenarios"
 
 * --- Sets used for the changing unit aggregation and efficiency approximations
     uft(unit, f, t) "Active units on intervals, enables aggregation of units for later intervals"

inc/1c_parameters.gms

@@ -126,6 +126,7 @@ Parameters
     df_central(f, t) "Displacement needed to reach the central forecast - this is needed when the forecast tree gets reduced in dynamic equations"
     df_reserves(node, restype, f, t) "Forecast index displacement needed to reach the realized forecast when committing reserves"
     df_scenario(f, t) "Forecast index displacement needed to get central forecast data for long-term scenarios"
+    df_realization(f, t) "Displacement needed to reach the realized forecast on the current time step when no forecast is available"
 
     // Sample displacement arrays
     ds(s, t) "Displacement needed to reach the sample of previous time step"

inc/1e_inputs.gms

@@ -581,3 +581,9 @@ loop(timeseries$(not sameas(timeseries, 'ts_cf')),
 p_autocorrelation(flowNode, 'ts_cf') = 0;
 p_tsMinValue(flowNode, 'ts_cf') = 0;
 p_tsMaxValue(flowNode, 'ts_cf') = 1;
+
+* By default all nodes use forecasts for all timeseries
+gn_forecasts(gn, timeseries) = yes;
+gn_forecasts(flowNode, timeseries) = yes;
+gn_forecasts(restype, node, 'ts_reserveDemand') = yes;
+

inc/3b_periodicLoop.gms

@@ -354,11 +354,6 @@ loop(s_scenario(s, scenario)$(ord(s) > 1 and ord(scenario) > 1),
              = (ord(scenario) - 1) * mSettings(mSolve, 'scenarioLength');
     );
 
-    loop(gn_scenarios(grid, node, param_gnBoundaryTypes),
-      dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
-          = (ord(scenario) - 1) * mSettings(mSolve, 'scenarioLength');
-    );
-
     loop(gn_scenarios(flow, node, timeseries),
         dt_scenarioOffset(flow, node, timeseries, s)
             = (ord(scenario) - 1) * mSettings(mSolve, 'scenarioLength');
@@ -448,6 +443,13 @@ df(f_solve(f), t_active(t))${ ord(t) <= tSolveFirst + max(mSettings(mSolve, 't_j
                                                           min(mSettings(mSolve, 't_perfectForesight'),
                                                               currentForecastLength))}
     = sum(mf_realization(mSolve, f_), ord(f_) - ord(f));
+// Displacement to reach the realized forecast
+loop(mf_realization(mSolve, f_),
+    df_realization(f_solve(f), t)$[ord(t) >= tSolveFirst
+                                   and ord(t) <= tSolveFirst + currentForecastLength
+                                  ]
+      = ord(f_) - ord(f);
+);
 // Central forecast for the long-term scenarios comes from a special forecast label
 if(mSettings(mSolve, 'scenarios') > 1,
     loop((ms_initial(mSolve, s), mf_central(mSolve, f)),
@@ -455,6 +457,14 @@ if(mSettings(mSolve, 'scenarios') > 1,
           = sum(mf_scenario(mSolve, f_), ord(f_) - ord(f));
     );
 );
+// Check that df_forecast and df_scenario do not overlap
+loop(ft(f, t),
+  if(df_realization(f, t) and df_scenario(f, t),
+      put log "!!! Overlapping period of using realization and scenarios"/;
+      put log "!!! Check forecast lengths, `gn_scenarios` and `gn_forecasts`"/;
+      abort "Overlapping realization and scenarios!";
+  );
+);
 
 // Forecast displacement between central and forecasted intervals at the end of forecast horizon
 Option clear = df_central; // This can be reset.

inc/3c_inputsLoop.gms

@@ -322,17 +322,19 @@ $offtext
     ts_influx_(gn(grid, node), fts(f, tt_interval(t), s))
         = sum(tt_aggregate(t, t_),
             ts_influx(grid, node,
-                f + (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'))))
+                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_cf_(flowNode(flow, node), fts(f, tt_interval(t), s))
         = sum(tt_aggregate(t, t_),
             ts_cf(flow, node,
-                f + (df_scenario(f, t)$gn_scenarios(flow, node, 'ts_cf')),
-                t_+ (   + dt_scenarioOffset(flow, node, 'ts_cf', s)
-                        + dt_circular(t_)$(not gn_scenarios(flow, node, 'ts_cf'))))
+                f + (  df_realization(f, t)$(not gn_forecasts(flow, node, 'ts_cf'))
+                     + df_scenario(f, t)$gn_scenarios(flow, node, 'ts_cf')),
+                t_+ (  dt_scenarioOffset(flow, node, 'ts_cf', s)
+                     + dt_circular(t_)$(not gn_scenarios(flow, node, 'ts_cf'))))
             )
             / mInterval(mSolve, 'stepsPerInterval', counter);
     // Reserves relevant only until reserve_length
@@ -340,7 +342,8 @@ $offtext
       ${ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')  }
         = sum(tt_aggregate(t, t_),
             ts_reserveDemand(restype, up_down, node,
-                f + (df_scenario(f, t)$gn_scenarios(restype, node, 'ts_reserveDemand')),
+                f + (  df_realization(f, t)$(not gn_forecasts(restype, node, 'ts_reserveDemand'))
+                     + df_scenario(f, t)$gn_scenarios(restype, node, 'ts_reserveDemand')),
                 t_+ dt_circular(t_))
             )
             / mInterval(mSolve, 'stepsPerInterval', counter);
@@ -349,7 +352,8 @@ $offtext
            // Take average if not a limit type
         = (sum(tt_aggregate(t, t_),
                 ts_node(grid, node, param_gnBoundaryTypes,
-                    f + (df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
+                    f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_node'))
+                         + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
                     t_+ (   + dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                             + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node'))))
             )
@@ -359,7 +363,8 @@ $offtext
           // Maximum lower limit
           + smax(tt_aggregate(t, t_),
                 ts_node(grid, node, param_gnBoundaryTypes,
-                    f + (df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
+                    f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_node'))
+                         + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
                     t_+ (   + dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                             + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node'))))
                 )
@@ -367,7 +372,8 @@ $offtext
           // Minimum upper limit
           + smin(tt_aggregate(t, t_),
                 ts_node(grid, node, param_gnBoundaryTypes,
-                    f + (df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
+                    f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_node'))
+                         + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
                     t_+ (   + dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                             + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node'))))
                 )

inc/smoothing.gms

@@ -31,7 +31,10 @@ $endif
     = min(p_tsMaxValue(%linking_set%, '%timeseries%'),
           max(p_tsMinValue(%linking_set%, '%timeseries%'),
               %timeseries%_(%linking_set%, f, t, s)
-              + (%timeseries%(%linking_set%, f_, t_)
+              + (%timeseries%(%linking_set%,
+                              f_ + (df_realization(f_, t_)
+                                    $(not gn_forecasts(%linking_set%, '%timeseries%'))),
+                              t_)
                  - %timeseries%(%linking_set%,
                              f + (df_scenario(f, t)$gn_scenarios(%linking_set%,
                                                                  '%timeseries%')),