Fixed the scenario smoothing algorithm

Standard deviations of scenario data now needs to be given as input.

inc/1c_parameters.gms

@@ -170,18 +170,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/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/3c_inputsLoop.gms

@@ -475,73 +475,17 @@ loop((mst_start(mSolve, s, t), ss(s, s_)),
 
 * --- Smooting of stochastic scenarios ----------------------------------------
 $ontext
-First calculate standard deviation for over all samples, then smoothen the scenarios
-following the methodology presented in [1, p. 443]. This avoids a discontinuity
-`jump' after the initial sample.
+Smoothen the scenarios following the methodology presented in [1, p. 443].
+This avoids a discontinuity `jump' after the initial sample.
 
 [1] A. Helseth, B. Mo, A. Lote Henden, and G. Warland, "Detailed long-term hydro-
     thermal scheduling for expansion planning in the Nordic power system," IET Gener.
     Transm. Distrib., vol. 12, no. 2, pp. 441 - 447, 2018.
 $offtext
 
-* Only do smoothing if scenarios are used
-* TODO: Also enable when long-term forecast is used
-if(mSettings(msolve, 'scenarios'),
-
-* Influx
-loop(gn(grid, node)$p_autocorrelation(grid, node, 'ts_influx'),
-    ts_influx_mean(grid, node, t_active(t))$(sum(sft(s, f, t), 1))
-        = sum(sft(s, f, t), ts_influx_(grid, node, f, t, s))
-                / sum(sft(s, f, t), 1);
-
-    ts_influx_std(grid, node, t_active(t))$(sum(sft(s, f, t), 1))
-        = sqrt(sum(sft(s, f, t), sqr(ts_influx_(grid, node, f, t, s)
-                                     - ts_influx_mean(grid, node, t)))
-                / sum(sft(s, f, t), 1)
-          );
-
-    // Do smoothing
-    loop(mst_end(ms_initial(mSolve, s_), t_),
-        ts_influx_(grid, node, ft(f, t), s)$(ts_influx_std(grid, node, t_+dt_circular(t_))
-                                             and sft(s, f, t)
-                                             and not ms_initial(mSolve, s))
-            = min(p_tsMaxValue(node, 'ts_influx'), max(p_tsMinValue(node, 'ts_influx'),
-              ts_influx_(grid, node, f, t, s)
-              + (ts_influx_(grid, node, f, t_, s_)
-                 - ts_influx_(grid, node, f, t_, s))
-                * (ts_influx_std(grid, node, t+dt_circular(t))
-                    / ts_influx_std(grid, node, t_+dt_circular(t_)))
-                * power(p_autocorrelation(grid, node, 'ts_influx'), abs(ord(t) - ord(t_)))
-              ));
-    );
+// Do smoothing
+loop((ms_initial(mSolve, s_), mf_central(mSolve, f_), t_)
+    $(ord(t_) = msEnd(mSolve, s_)),
+    $$batinclude 'inc/smoothing.gms' ts_influx
+    $$batinclude 'inc/smoothing.gms' ts_cf
 );
-
-* CF
-loop(flowNode(flow, node)$p_autocorrelation(flow, node, 'ts_cf'),
-    ts_cf_mean(flow, node, t_active(t))$(sum(sft(s, f, t), 1))
-        = sum(sft(s, f, t), ts_cf_(flow, node, f, t, s))
-                / sum(sft(s, f, t), 1);
-
-    ts_cf_std(flow, node, t_active(t))$(sum(sft(s, f, t), 1))
-        = sqrt(sum(sft(s, f, t), sqr(ts_cf_(flow, node, f, t, s)
-                                     - ts_cf_mean(flow, node, t)))
-                / sum(sft(s, f, t), 1)
-          );
-
-    // Do smoothing
-    loop((mst_end(ms_initial(mSolve, s_), t_), f_)$ft(f_, t_),
-        ts_cf_(flow, node, ft(f, t), s)$(ts_cf_std(flow, node, t_+dt_circular(t_))
-                                         and sft(s, f, t)
-                                         and not ms_initial(mSolve, s))
-            = min(p_tsMaxValue(node, 'ts_cf'), max(p_tsMinValue(node, 'ts_cf'),
-              ts_cf_(flow, node, f, t, s)
-              + (ts_cf_(flow, node, f_, t_, s_)
-                 - ts_cf_(flow, node, f, t_, s))
-                * (ts_cf_std(flow, node, t+dt_circular(t))
-                    / ts_cf_std(flow, node, t_+dt_circular(t_)))
-                * power(p_autocorrelation(flow, node, 'ts_cf'), abs(ord(t) - ord(t_)))
-              ));
-    );
-);
-
-); // end if(mSettings(msolve, 'scenarios'),

inc/smoothing.gms

+$ontext
+This file is part of Backbone.
+
+Backbone is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+Backbone is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+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
+
+* Argument gives time series type
+$setargs timeseries
+
+* Select set linking nodes based on the time series type
+$ifthen %timeseries% == 'ts_cf'
+$setlocal linking_set flowNode
+$else
+$setlocal linking_set gn
+$endif
+
+%timeseries%_(%linking_set%, fts(f, t, s))$[p_autocorrelation(%linking_set%,
+                                                              '%timeseries%')
+                                            and ms_central(mSolve, s)]
+    = 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_scenario(f, t)$gn_scenarios(%linking_set%,
+                                                                 '%timeseries%')),
+                             t_ + dt_scenarioOffset(%linking_set%,
+                                                    '%timeseries%', s)))
+                * [(%timeseries%_std(%linking_set%, t+dt_circular(t))
+                    / %timeseries%_std(%linking_set%, t_+dt_circular(t_))
+                   )$%timeseries%_std(%linking_set%, t_+dt_circular(t_))
+                   // If standard deviations not defined, use 1
+                   + 1$(not %timeseries%_std(%linking_set%, t_+dt_circular(t_)))
+                  ]
+                * power(p_autocorrelation(%linking_set%, '%timeseries%'),
+                        abs(ord(t) - ord(t_)))
+      ));