Merge branch 'dev' into nhniina_rampCosts

Resolved conflicts in inc/2d_constraints.gms (e.g. start-up and shutdown variables).

4d_postProcess_invest2schedule_temp.gms

+* =============================================================================
+* --- Save investments results in changes.inc file to be used in child setups -
+* =============================================================================
+
+* Output file streams
+file f_changes /'output\changes.inc'/;
+
+put f_changes
+
+* Do not allow investments in the child setups
+loop(unit,
+    put "p_unit('", unit.tl, "', 'maxUnitCount') = 0;"/;
+);
+
+* Update the number of subunits in the child setups (rounded here to the nearest integer)
+loop(unit${r_invest(unit)},
+    tmp = round(r_invest(unit), 0)
+    put "p_unit('", unit.tl, "', 'unitCount') = p_unit('", unit.tl, "', 'unitCount') + ", tmp, ";"/;
+);
+* Update maxGen and maxCons values in the child setups
+loop(gnu(grid, node, unit)${r_invest(unit)},
+    tmp = round(r_invest(unit), 0) * p_gnu(grid, node, unit, 'unitSizeGen');
+    put "p_gnu('", grid.tl, "', '", node.tl, "', '", unit.tl, "', 'maxGen') = p_gnu('", grid.tl, "', '", node.tl, "', '", unit.tl, "', 'maxGen') + ", tmp, ";"/;
+    tmp = round(r_invest(unit), 0) * p_gnu(grid, node, unit, 'unitSizeCons');
+    put "p_gnu('", grid.tl, "', '", node.tl, "', '", unit.tl, "', 'maxCons') = p_gnu('", grid.tl, "', '", node.tl, "', '", unit.tl, "', 'maxCons') + ", tmp, ";"/;
+);
+
+* Do not allow investments in the child setups (commented out at the moment)
+*loop(gn2n_directional(grid, node, node_),
+*    put "p_gnn('", grid.tl, "', '", node.tl, "', '", node_.tl, "', 'transferCapMax') = 0;"/;
+*    put "p_gnn('", grid.tl, "', '", node_.tl, "', '", node.tl, "', 'transferCapMax') = 0;"/;
+*);
+
+* Update transmission capacity in the child setups
+loop(gn2n_directional(grid, node, node_)${sum(t_invest, r_investTransfer(grid, node, node_, t_invest))},
+    tmp = sum(t_invest, r_investTransfer(grid, node, node_, t_invest));
+    put "p_gnn('", grid.tl, "', '", node.tl, "', '", node_.tl, "', 'transferCap') = p_gnn('", grid.tl, "', '", node.tl, "', '", node_.tl, "', 'transferCap') + ", tmp, ";"/;
+    put "p_gnn('", grid.tl, "', '", node_.tl, "', '", node.tl, "', 'transferCap') = p_gnn('", grid.tl, "', '", node_.tl, "', '", node.tl, "', 'transferCap') + ", tmp, ";"/;
+);
+
+putclose;

CHANGELOG.md

@@ -2,6 +2,33 @@
 All notable changes to this project will be documented in this file.
 
 ## [Unreleased]
+## Added
+- New model setting 't_perfectForesight' tells the number of time steps (from 
+  the beginning of current solve) for which realized data is used instead of 
+  forecasts. This value cannot exceed current forecast length, however. Setting 
+  the value lower than 't_jump' has no effect.
+- Automated the calculation of sample start and end times if using long-term 
+  scenarios. Also setting number of scenarios to one, instructs the model to use
+  central forecast for the long-term.
+- Speedup for model dimension calculation (set `msft` etc.)
+
+## [1.0.6] - 2019-03-27
+### Fixed
+- Major bug in state variable reserve equations
+- Scenario smoothing alogirithm
+
+### Changed
+- Speedup for timeseries calculations
+
+### Added 
+- New model setting `mSettings(mType, 'onlyExistingForecasts') = 0|1` to control 
+  the reading of forecasts. Set to 1 to only read forecast data that exists in 
+  the file. Note that zeros need to be saved as Eps when using this.
+- Proper stochastic programming for the long-term scenarios period. Possible also
+  to create a stochastic tree from the original data.
+- Clickable link to *sr.log* in the process window in case of SCENRED2 error
+- New diagnostic parameter for timeseries scenarios `d_ts_scenarios`
+
 
 ## [1.0.5] - 2019-02-14
 ### Fixed
@@ -55,7 +82,8 @@ All notable changes to this project will be documented in this file.
 
 
 
-[Unreleased]: https://gitlab.vtt.fi/backbone/backbone/compare/v1.0.5...dev
+[Unreleased]: https://gitlab.vtt.fi/backbone/backbone/compare/v1.0.6...dev
+[1.0.6]: https://gitlab.vtt.fi/backbone/backbone/compare/v1.0.5...v1.0.6
 [1.0.5]: https://gitlab.vtt.fi/backbone/backbone/compare/v1.0.4...v1.0.5
 [1.0.4]: https://gitlab.vtt.fi/backbone/backbone/compare/v1.0.3...v1.0.4
 [1.0.3]: https://gitlab.vtt.fi/backbone/backbone/compare/v1.0.2...v1.0.3

defModels/investInit_temp.gms

@@ -41,11 +41,9 @@ if (mType('invest'),
     // Number of samples used by the model
     mSettings('invest', 'samples') = 3;
 
-    // Define Initial and Central samples
+    // Clear Initial and Central samples
     ms_initial('invest', s) = no;
-    ms_initial('invest', 's000') = yes;
     ms_central('invest', s) = no;
-    ms_central('invest', 's000') = yes;
 
     // Define time span of samples
     msStart('invest', 's000') = 1;

defModels/scheduleInit_temp.gms

@@ -39,47 +39,45 @@ if (mType('schedule'),
 * --- Model Time Structure ----------------------------------------------------
 * =============================================================================
 
+* --- Long-term scenarios -----------------------------------------------------
+
+    // Number of long-term scenarios to use (zero if none)
+    mSettings('schedule', 'scenarios') = 0;
+
+    // Length of single scenario (in time steps)
+    mSettings('schedule', 'scenarioLength') = 8760;
+
 * --- Define Samples ----------------------------------------------------------
 
     // Number of samples used by the model
     mSettings('schedule', 'samples') = 1;
 
-    // Sample length of stocahstic data
-    mSettings('schedule', 'sampleLength') = mSettings('schedule', 'dataLength');
-
     // Define Initial and Central samples
     ms_initial('schedule', s) = no;
-    ms_initial('schedule', 's000') = yes;
-    ms_central('schedule', s) = no;
-    ms_central('schedule', 's000') = yes;
+    ms_central('schedule', s) = no; 
 
     // Define time span of samples
     msStart('schedule', 's000') = 1;
     msEnd('schedule', 's000') = msStart('schedule', 's000') + mSettings('schedule', 't_horizon');
 
-    // If using long-term samples, uncomment
-    //msEnd('schedule', 's000') = msStart('schedule', 's000')
-    //                            + mSettings('schedule', 't_forecastLengthUnchanging');
-
-    //loop(s $(ord(s) > 1),
-    //    msStart('schedule', s) = msStart('schedule', 's000')
-    //                             + mSettings('schedule', 't_forecastLengthUnchanging');
-    //    msEnd('schedule', s) = msStart('schedule', 's000')
-    //                           + mSettings('schedule', 't_horizon');
-    //);
-
     // Define the probability (weight) of samples
     p_msProbability('schedule', s) = 0;
     p_msProbability('schedule', 's000') = 1;
 
+    // If using long-term samples, uncomment
+    //ms_central('schedule', 's001') = yes;
+    //
+    //msEnd('schedule', 's000') = msStart('schedule', 's000') + 168;
+
+    //msStart('schedule', 's001') = msEnd('schedule', 's000');
+    //msEnd('schedule', 's001') = msStart('schedule', 's000')
+    //                       + mSettings('schedule', 't_horizon');
+    //
+    //p_msProbability('schedule', 's001') = 1;
+    //);
+
     // Define which nodes use long-term samples
-    loop(gn(grid, node)$sameas(grid, 'hydro'),
-        //longtermSamples('hydro', node, 'ts_influx') = yes;
-        //longtermSamples('hydro', node, 'minSpill') = yes;
-    );
-    loop(flowNode(flow, node)$sameas(flow, 'wind'),
-        //longtermSamples('wind', node, 'ts_cf') = yes;
-    );
+    //gn_scenarios('hydro', 'XXX', 'ts_influx') = yes;
 
 * --- Define Time Step Intervals ----------------------------------------------
 
@@ -109,6 +107,7 @@ if (mType('schedule'),
     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)
     mSettings('schedule', 't_forecastLengthDecreasesFrom') = 168;  // Length of forecasts in time steps - this decreases when the solve moves forward until the new forecast data is read (then extends back to full length)
+    mSettings('schedule', 't_perfectForesight') = 0;               // How many time steps after there is perfect foresight (including t_jump)
     mSettings('schedule', 't_forecastJump') = 24;                  // How many time steps before new forecast is available
     mSettings('schedule', 't_improveForecast') = 0;                // Number of time steps ahead of time that the forecast is improved on each solve.
 

defOutput/debugSymbols.inc

@@ -79,8 +79,10 @@ v_sos2
 v_spill
 v_online_LP
 v_online_MIP
-v_startup
-v_shutdown
+v_startup_LP
+v_startup_MIP
+v_shutdown_LP
+v_shutdown_MIP
 v_genRampUpDown
 v_transferRightward
 v_transferLeftward

inc/1a_definitions.gms

@@ -50,6 +50,8 @@ Sets
         t_forecastStart, // Time step for first reading the forecasts (not necessarily t_start)
         t_forecastJump, // Number of time steps between each update of the forecasts
         t_improveForecast "Number of time steps ahead of time on which the forecast is improved on each solve"
+        t_perfectForesight "Number of time steps for which realized data is used instead of forecasts"
+        onlyExistingForecasts "Use only existing forecast values when reading updated forecasts. Note: zero values need to be saved as Eps in the gdx file."
         scenarios        "Number of long-term scenarios used"
         scenarioLength   "Length of scenario in time steps for creating stocahstic scenarios from time series data"
 
@@ -378,9 +380,9 @@ hr(param_unit) "Heat rate for the corresponding operating point ('hrop') in the
 // Slack categories
 slack(param_gnBoundaryTypes) "Categories for slack variables"
        / upwardSlack01*upwardSlack20, downwardSlack01*downwardSlack20 /
-upwardSlack(slack) "Set of upward slacks"
+upwardSlack(param_gnBoundaryTypes) "Set of upward slacks"
        / upwardSlack01*upwardSlack20 /
-downwardSlack(slack) "Set of downward slacks"
+downwardSlack(param_gnBoundaryTypes) "Set of downward slacks"
        / downwardSlack01*downwardSlack20 /
 
 // Flags for boundaries

inc/1b_sets.gms

@@ -94,6 +94,7 @@ Sets
     tt_block(counter, t) "Temporary time step subset for storing the time interval blocks"
     tt_interval(t) "Temporary time steps when forming the ft structure, current sample"
     tt_forecast(t) "Temporary subset for time steps used for forecast updating during solve loop"
+    tt_aggregate(t, t) "Time steps included in each active time step for time series aggregation"
     mf(mType, f) "Forecasts present in the models"
     ms(mType, s) "Samples present in the models"
     mst_start(mType, s, t) "Start point of samples"
@@ -106,12 +107,14 @@ Sets
     msf(mType, s, f) "Combination of samples and forecasts in the models"
     msft(mType, s, f, t) "Combination of models, samples, forecasts and t's"
     sft(s, f, t) "Combination of samples, forecasts and t's in the current model solve"
+    fts(f, t, s) "Like `sft` but different order"
     sft_realized(s, f, t) "Realized sft"
     sft_realizedNoReset(s, f, t) "Full set of realized sft, facilitates calculation of results"
     msft_realizedNoReset(mType, s, f, t) "Combination of realized samples, forecasts and t:s in the current model solve and previously realized t:s"
     mft_start(mType, f, t) "Start point of the current model solve"
-    mf_realization(mType, f) "fRealization of the forecasts"
+    mf_realization(mType, f) "Realization of the forecasts"
     mf_central(mType, f) "Forecast that continues as sample(s) after the forecast horizon ends"
+    mf_scenario(mType, f) "Forecast label that is used for long-term scenario data"
     ms_initial(mType, s) "Sample that presents the realized/forecasted period"
     ms_central(mType, s) "Sample that continues the central forecast after the forecast horizon ends"
     mft_lastSteps(mType, f, t) "Last interval of the current model solve"
@@ -121,9 +124,9 @@ Sets
     gnss_bound(grid, node, s, s) "Bound the samples so that the node state at the last interval of the first sample equals the state at the first interval of the second sample"
     uss_bound(unit, s, s) "Bound the samples so that the unit online state at the last interval of the first sample equals the state at the first interval of the second sample"
     s_active(s) "Samples with non-zero probability in the current model solve"
-    s_stochastic(s) "Samples used for stochastic programming"
     ss(s, s) "Previous sample of sample"
     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"
 
@@ -167,6 +170,7 @@ Option clear = modelSolves;
 Option clear = ms;
 Option clear = mf;
 mf_realization(mType, 'f00') = yes;
+Option clear = mf_scenario;
 Option clear = gn_scenarios;
 Option clear = mTimeseries_loop_read;
 

inc/1c_parameters.gms

@@ -97,9 +97,13 @@ Parameters
     p_mfProbability(mType, f) "Probability of forecast"
     p_msft_probability(mType, s, f, t) "Probability of forecast"
     p_sProbability(s) "Probability of sample"
+$if defined scenario
     p_scenProbability(scenario) "Probability of scenarios"
 ;
 
+$if declared p_scenProbability
+Option clear = p_scenProbability;  // Initialize with empty data
+
 Scalar p_sWeightSum "Sum of sample weights";
 
 * --- Model structure ---------------------------------------------------------
@@ -108,7 +112,7 @@ Parameters
     dt(t) "Displacement needed to reach the previous time interval (in time steps)"
     dt_circular(t) "Circular t displacement if the time series data is not long enough to cover the model horizon"
     dt_next(t) "Displacement needed to reach the next time interval (in time steps)"
-    dtt(t, t) "Displacement needed to reach any previous time interval (in time steps)"
+    dt_active(t) "Displacement needed to reach the corresponding active time interval from any time interval (in time steps)"
     dt_toStartup(unit, t) "Displacement from the current time interval to the time interval where the unit was started up in case online variable changes from 0 to 1 (in time steps)"
     dt_toShutdown(unit, t) "Displacement from the current time interval to the time interval where the shutdown phase began in case generation becomes 0 (in time steps)"
     dt_starttypeUnitCounter(starttype, unit, counter) "Displacement needed to account for starttype constraints (in time steps)"
@@ -121,7 +125,7 @@ Parameters
     df(f, t) "Displacement needed to reach the realized forecast on the current time step"
     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 realized scenario data for long-term scenarios"
+    df_scenario(f, t) "Forecast index displacement needed to get central forecast data for long-term scenarios"
 
     // Sample displacement arrays
     ds(s, t) "Displacement needed to reach the sample of previous time step"
@@ -133,7 +137,7 @@ Parameters
     ddf_(f) "Temporary forecast displacement array"
 
     // Other
-    p_slackDirection(slack) "+1 for upward slacks and -1 for downward slacks"
+    p_slackDirection(param_gnBoundaryTypes) "+1 for upward slacks and -1 for downward slacks"
     tForecastNext(mType) "When the next forecast will be available (ord time)"
     aaSolveInfo(mType, t, solveInfoAttributes) "Stores information about the solve status"
     msStart(mType, s) "Start point of samples: first time step in the sample"
@@ -172,18 +176,19 @@ Parameters
     ts_unavailability_update(unit, t)
 
     // Help parameters for calculating smoothening of time series
-    ts_influx_mean(grid, node, t) "Mean of ts_influx over samples"
     ts_influx_std(grid, node, t)  "Standard deviation of ts_influx over samples"
-    ts_cf_mean(flow, node, t) "Mean of ts_cf over samples (p.u.)"
     ts_cf_std(flow, node, t) "Standard deviation of ts_cf over samples (p.u.)"
 
     p_autocorrelation(*, node, timeseries) "Autocorrelation of time series for the grid/flow, node and time series type (lag = 1 time step)"
 
     // Bounds for scenario smoothening
-    p_tsMinValue(node, timeseries)
-    p_tsMaxValue(node, timeseries)
+    p_tsMinValue(*, node, timeseries) "Minimum allowed value of timeseries for grid/flow and node"
+    p_tsMaxValue(*, node, timeseries) "Maximum allowed value of timeseries in grid/flow and node"
 ;
 
+* Reset values for some parameters
+Options clear = ts_influx_std, clear = ts_cf_std;
+
 * --- Other time dependent parameters -----------------------------------------
 Parameters
     p_storageValue(grid, node, t) "Value of stored something at the end of a time step"

inc/1d_results.gms

@@ -194,5 +194,6 @@ Parameters
     d_nodeState(grid, node, param_gnBoundaryTypes, s, f, t) "Diagnostic temperature forecasts (accounting for GAMS plotting error)"
     d_influx(grid, node, s, f, t) "Diagnostic influx forecasts (accounting for GAMS plotting error)"
     d_state(grid, node, scenario, f, t) "Diagnostic state results in each scenario"
+    d_ts_scenarios(timeseries, *, node, scenario, f, t) "Diagnostic time series values in scenarios"
 ;
 $endif.diag

inc/1e_inputs.gms

@@ -535,9 +535,11 @@ loop( (nu(node, unit), restypeDirection(restype, up_down)),
 * =============================================================================
 * --- Default values  ---------------------------------------------------------
 * =============================================================================
-
-loop((gn(grid, node), timeseries),
-    p_autocorrelation(grid, node, timeseries) = 0;
-    p_tsMinValue(node, timeseries) = -Inf;
-    p_tsMaxValue(node, timeseries) = Inf;
+loop(timeseries$(not sameas(timeseries, 'ts_cf')),
+    p_autocorrelation(gn, timeseries) = 0;
+    p_tsMinValue(gn, timeseries) = -Inf;
+    p_tsMaxValue(gn, timeseries) = Inf;
 );
+p_autocorrelation(flowNode, 'ts_cf') = 0;
+p_tsMinValue(flowNode, 'ts_cf') = 0;
+p_tsMaxValue(flowNode, 'ts_cf') = 1;

inc/2a_variables.gms

@@ -23,6 +23,8 @@ Free variables
     v_transfer(grid, node, node, s, f, t) "Average electricity transmission level from node to node during an interval (MW)"
 ;
 Integer variables
+    v_startup_MIP(unit, starttype, s, f, t) "Sub-units started up after/during an interval (p.u.), (MIP variant)"
+    v_shutdown_MIP(unit, s, f, t) "Sub-units shut down after/during an interval (p.u.) (MIP variant)"
     v_online_MIP(unit, s, f, t) "Number of sub-units online for units with unit commitment restrictions"
     v_invest_MIP(unit, t) "Number of invested sub-units"
     v_investTransfer_MIP(grid, node, node, t) "Number of invested transfer links"
@@ -35,8 +37,8 @@ SOS2 variables
 ;
 Positive variables
     v_fuelUse(fuel, unit, s, f, t) "Fuel use of a unit during an interval (MWh_fuel)"
-    v_startup(unit, starttype, s, f, t) "Sub-units started up after/during an interval (p.u.)"
-    v_shutdown(unit, s, f, t) "Sub-units shut down after/during an interval (p.u.)"
+    v_startup_LP(unit, starttype, s, f, t) "Sub-units started up after/during an interval (p.u.), (LP variant)"
+    v_shutdown_LP(unit, s, f, t) "Sub-units shut down after/during an interval (p.u.) (LP variant)"
     v_genRampUpDown(grid, node, unit, slack, s, f, t) "Change in energy generation or consumption over an interval, separated into different 'slacks' (MW/h)"
     v_spill(grid, node, s, f, t) "Spill of energy from storage node during an interval (MWh)"
     v_transferRightward(grid, node, node, s, f, t) "Average electricity transmission level from the first node to the second node during an interval (MW)"

inc/2c_objective.gms

@@ -89,7 +89,10 @@ q_obj ..
                 // Start-up costs, initial startup free as units could have been online before model started
                 + sum(uft_online(unit, f, t),
                     + sum(unitStarttype(unit, starttype),
-                        + v_startup(unit, starttype, s, f, t) // Cost of starting up
+                        + [ // Unit startup variables
+                            + v_startup_LP(unit, starttype, s, f, t)${ uft_onlineLP(unit, f, t) }
+                            + v_startup_MIP(unit, starttype, s, f, t)${ uft_onlineMIP(unit, f, t) }
+                            ]
                             * [ // Startup variable costs
                                 + p_uStartup(unit, starttype, 'cost')
 
@@ -111,7 +114,13 @@ q_obj ..
 
                 // Shut-down costs, initial shutdown free?
                 + sum(uft_online(unit, f, t),
-                      + v_shutdown(unit, s, f, t) * p_uShutdown(unit, 'cost')
+                    + p_uShutdown(unit, 'cost')
+                      * [
+                            + v_shutdown_LP(unit, s, f, t)
+                                ${ uft_onlineLP(unit, f, t) }
+                            + v_shutdown_MIP(unit, s, f, t)
+                                ${ uft_onlineMIP(unit, f, t) }
+                        ]
                   ) // END sum(uft_online)
 
                 // Ramping costs

inc/3a_periodicInit.gms

@@ -53,11 +53,23 @@ $offtext
     // Set the time for the next available forecast.
     tForecastNext(m) = mSettings(m, 't_forecastStart');
 
+    // Update number of samples
+    if(mSettings(m, 'scenarios') = 1, mSettings(m, 'samples') = 2);
+
     // Select samples for the model
     if (not sum(s, ms(m, s)),  // unless they have been provided as input
         ms(m, s)$(ord(s) <= mSettings(m, 'samples')) = yes;
     );
 
+    // Set active and previous samples
+    loop(ms(m, s),
+        s_active(s) = yes;
+        loop(s_ $ms(m, s_),
+            // Set previous samples for samples
+            ss(s, s_)$(msStart(m, s) = msEnd(m, s_)) = yes;
+        );
+    );
+
     // Select forecasts in use for the models
     if (not sum(f, mf(m, f)),  // unless they have been provided as input
         mf(m, f)$(ord(f) <= 1 + mSettings(m, 'forecasts')) = yes;  // realization needs one f, therefore 1 + number of forecasts
@@ -69,6 +81,9 @@ $offtext
 
     // Select combinations of models, samples and forecasts to be solved
     msf(m, s, f_solve(f))$(ms(m, s) and mf(m, f)) = yes;
+    if(mSettings(m, 'scenarios'),
+        msf(ms_central(m, s), f_solve(f))$mf_realization(m, f) = no;
+    );
 
     // Check the modelSolves for preset patterns for model solve timings
     // If not found, then use mSettings to set the model solve timings
@@ -628,6 +643,15 @@ loop(msf(m, s, f)${ mf_realization(m, f) },
 
 loop(m, // Not ideal, but multi-model functionality is not yet implemented
 
+* --- Prefect foresight not longer than forecast length
+    if(mSettings(m, 't_perfectForesight')
+       > max(mSettings(m, 't_forecastLengthUnchanging'),
+             mSettings(m, 't_forecastLengthDecreasesFrom')),
+        put log "!!! Error in model ", m.tl:0 /;
+        put log "!!! Error: t_perfectForesight > max(t_forecastLengthUnchanging, t_forecastLengthDecreasesFrom)"/;
+        abort "Period of perfect foresight cannot be longer than forecast horizon";
+    );
+
 * --- Reserve structure checks ------------------------------------------------
 
     loop(restypeDirectionNode(restype, up_down, node),

inc/3b_periodicLoop.gms

@@ -36,8 +36,10 @@ Option clear = v_help_inc;
 Option clear = v_sos2;
 // Positive Variables
 Option clear = v_fuelUse;
-Option clear = v_startup;
-Option clear = v_shutdown;
+Option clear = v_startup_LP;
+Option clear = v_startup_MIP;
+Option clear = v_shutdown_LP;
+Option clear = v_shutdown_MIP;
 Option clear = v_genRampUpDown;
 Option clear = v_spill;
 Option clear = v_transferRightward;
@@ -145,12 +147,19 @@ Option clear = p_stepLength;
 Option clear = msft;
 Option clear = mft;
 Option clear = ft;
-Option clear = sft;
-Option clear = s_active, clear = s_stochastic, clear = s_scenario, clear = ss;
+Option clear = sft, clear = fts;
 Option clear = mst_start, clear = mst_end;
+$ifthen declared scenario
+if(mSettings(mSolve, 'scenarios'),  // Only clear these if using long-term scenarios
+    Options clear = s_active, clear = s_scenario, clear = ss,
+            clear = p_msProbability, clear = ms_central;
+);
+$endif
+
 
 // Initialize the set of active t:s, counters and interval time steps
 Option clear = t_active;
+Option clear = dt_active;
 Option clear = tt_block;
 Option clear = cc;
 tCounter = 1;
@@ -207,14 +216,17 @@ loop(cc(counter),
 
     // If stepsPerInterval equals one, simply use all the steps within the block
     if(mInterval(mSolve, 'stepsPerInterval', counter) = 1,
-        tt_interval(tt(t)) = yes; // Include all time steps within the block
+        // Include all time steps within the block
+        tt_interval(tt(t)) = yes;
 
     // If stepsPerInterval exceeds 1 (stepsPerInterval < 1 not defined)
     elseif mInterval(mSolve, 'stepsPerInterval', counter) > 1,
-        tt_interval(tt(t)) // Select the active time steps within the block
-             ${mod(ord(t) - tSolveFirst - tCounter,
-                   mInterval(mSolve, 'stepsPerInterval', counter)) = 0
-              } = yes;
+
+        // Calculate the displacement required to reach the corresponding active time step from any time step
+        dt_active(tt(t)) = - (mod(ord(t) - tSolveFirst - tCounter, mInterval(mSolve, 'stepsPerInterval', counter)));
+
+        // Select the active time steps within the block
+        tt_interval(tt(t))${ not dt_active(t) } = yes;
 
     ); // END ELSEIF intervalLenght
 
@@ -225,76 +237,102 @@ loop(cc(counter),
 
     // Determine the combinations of forecasts and intervals
     // Include the t_jump for the realization
-    mft(mf(mSolve, f_solve), tt_interval(t))
-       ${ord(t) <= tSolveFirst + mSettings(mSolve, 't_jump')
+    ft(f_solve, tt_interval(t))
+       ${ord(t) <= tSolveFirst + max(mSettings(mSolve, 't_jump'),
+                                     min(mSettings(mSolve, 't_perfectForesight'),
+                                         currentForecastLength))
          and mf_realization(mSolve, f_solve)
         } = yes;
 
     // Include the full horizon for the central forecast
-    mft(mf(mSolve, f_solve), tt_interval(t))
-      ${ord(t) > tSolveFirst + mSettings(mSolve, 't_jump')
+    ft(f_solve, tt_interval(t))
+      ${ord(t) > tSolveFirst + max(mSettings(mSolve, 't_jump'),
+                                   min(mSettings(mSolve, 't_perfectForesight'),
+                                       currentForecastLength))
         and (mf_central(mSolve, f_solve)
-             or mSettings('schedule', 'forecasts') = 0)
+             or mSettings(mSolve, 'forecasts') = 0)
        } = yes;
 
     // Include up to forecastLength for remaining forecasts
-    mft(mf(mSolve, f_solve), tt_interval(t))
+    ft(f_solve, tt_interval(t))
       ${not mf_central(mSolve, f_solve)
         and not mf_realization(mSolve, f_solve)
-        and ord(t) > tSolveFirst + mSettings(mSolve, 't_jump')
+        and ord(t) > tSolveFirst + max(mSettings(mSolve, 't_jump'),
+                                       min(mSettings(mSolve, 't_perfectForesight'),
+                                           currentForecastLength))
         and ord(t) <= tSolveFirst + currentForecastLength
        } = yes;
 
-    // Set of locked combinations of forecasts and intervals for the reserves?
-
     // Update tActive
     t_active(tt_interval) = yes;
 
-    // Loop over defined samples
-    msft(msf(mSolve, s, f_solve), tt_interval(t))
-        ${ord(t) > msStart(mSolve, s) + tSolveFirst - 1 // Move the samples along with the dispatch
-          and ord(t) < msEnd(mSolve, s) + tSolveFirst   // Move the samples along with the dispatch
-         } = mft(mSolve, f_solve, t);
-
-    // Create stochastic programming scenarios
-    if(mSettings(mSolve, 'scenarios'),
-        // Select root sample and central forecast
-        loop((ms_initial(mSolve, s_), mf_central(mSolve, f_)),
-            s_active(s_) = yes;
-            loop(scenario$(ord(scenario) <= mSettings(mSolve, 'scenarios')),
-                s_scenario(s_, scenario) = yes;
-                loop(tt_interval(t)$(ord(t) >= msEnd(mSolve, s_) + tSolveFirst),
-                    mft(mSolve, f_, t) = yes;
-                    loop(s$(ord(s) = mSettings(mSolve, 'samples') + count_sample),
-                        s_active(s) = yes;
-                        s_stochastic(s) = yes;
-                        msft(mSolve, s, f_, t) = yes;
-                        s_scenario(s, scenario) = yes;
-                        p_msProbability(mSolve, s) = p_scenProbability(scenario);
-                        msStart(mSolve, s) = ord(t) - tSolveFirst;
-                        msEnd(mSolve, s)
-                          = ord(t) - tSolveFirst
-                            + mInterval(mSolve, 'stepsPerInterval', counter);
-                    );
-                    count_sample = count_sample + 1;
+    // Update tCounter for the next block of intervals
+    tCounter = mInterval(mSolve, 'lastStepInIntervalBlock', counter) + 1;
+
+); // END loop(counter)
+
+// Reset initial sample start and end times if using scenarios
+if(mSettings(mSolve, 'scenarios'),
+    Option clear = msStart, clear = msEnd;
+    msStart(ms_initial) = 1;
+    msEnd(ms_initial) = currentForecastLength + 1;
+);
+
+$ifthen defined scenario
+// Create stochastic programming scenarios
+// Select root sample and central forecast
+loop((ms_initial(mSolve, s_), mf_central(mSolve, f)),
+    s_active(s_) = yes;
+    p_msProbability(mSolve, s_) = 1;
+    loop(scenario$(ord(scenario) <= mSettings(mSolve, 'scenarios')),
+        s_scenario(s_, scenario) = yes;
+        if(mSettings(mSolve, 'scenarios') > 1,
+            loop(ft(f, t)$(ord(t) >= msEnd(mSolve, s_) + tSolveFirst),
+                loop(s$(ord(s) = mSettings(mSolve, 'samples') + count_sample),
+                    s_active(s) = yes;
+                    ms_central(mSolve, s) = yes;
+                    s_scenario(s, scenario) = yes;
+                    p_msProbability(mSolve, s) = p_scenProbability(scenario);
+                    msStart(mSolve, s) = ord(t) - tSolveFirst;
+                    msEnd(mSolve, s) = ord(t) - tSolveFirst
+                                              + p_stepLength(mSolve, f, t);
                 );
+                count_sample = count_sample + 1;
+            );
+        elseif mSettings(mSolve, 'scenarios') = 1,
+            loop(ms(mSolve, s)$(not sameas(s, s_)),
+                s_active(s) = yes;
+                ms_central(mSolve, s) = yes;
+                p_msProbability(mSolve, s) = 1;
+                s_scenario(s, scenario) = yes;
+                msStart(mSolve, s) = msEnd(mSolve, s_);
+                msEnd(mSolve, s) = msStart(mSolve, s_)
+                                   + mSettings(mSolve, 't_horizon');
             );
-            msf(mSolve, s, f_) = s_active(s);
         );
-        ms(mSolve, s) = s_active(s);
     );
-    // Reduce the model dimension
-    ft(f_solve, tt_interval(t)) = mft(mSolve, f_solve, t);
-
-    // Reduce the sample dimension
-    sft(s, f, t)$msft(mSolve, s, f, t) = ft(f, t);
-
-    // Update tCounter for the next block of intervals
-    tCounter = mInterval(mSolve, 'lastStepInIntervalBlock', counter) + 1;
+    ms(ms_central(mSolve, s)) = yes;
+    msf(ms_central(mSolve, s), f) = yes;
+);
+$endif
+
+// Loop over defined samples
+loop(msf(mSolve, s, f)$msStart(mSolve, s),
+                      // Move the samples along with the dispatch
+    sft(s, ft(f, t))${ord(t) > msStart(mSolve, s) + tSolveFirst - 1
+                      and ord(t) < msEnd(mSolve, s) + tSolveFirst
+                     } = yes;
+);
 
-); // END loop(counter)
+// Update the model specific sets and the reversed dimension set
+mft(mSolve, ft(f, t)) = yes;
+ms(mSolve, s)$ms(mSolve, s) = s_active(s);
+msf(mSolve, s, f)$msf(mSolve, s, f) = s_active(s);
+msft(mSolve, sft(s, f, t)) = yes;
+fts(ft(f, t), s)$sft(s, f, t) = yes;
 
 * Build stochastic tree by definfing previous samples
+$ifthen defined scenario