3c_inputsLoop.gms 22.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
$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
17

18
* =============================================================================
19
* --- Update the Forecast Data ------------------------------------------------
20
* =============================================================================
21

22
put log 'ord tSolve: ';
23
put log ord(tSolve):0:0 /;
24
putclose log;
25

26
27
28
29
30
31
32
33
34
35
36
// Determine the necessary horizon for updating data
option clear = tmp;
tmp = max(  mSettings(mSolve, 't_forecastLengthUnchanging') + mSettings(mSolve, 't_forecastJump'),
            mSettings('schedule', 't_forecastLengthDecreasesFrom')
            );

// Find time steps until the forecast horizon
option clear = tt_forecast;
tt_forecast(t_current(t))
    ${ ord(t) <= tSolveFirst + tmp }
    = yes;
37

38
if (ord(tSolve) >= tForecastNext(mSolve),
39
$ontext
40
* These don't work due to the wild cards in the parameter definitions!
41
42
43
    // Update ts_unit
    if (mTimeseries_loop_read(mSolve, 'ts_unit'),
        put_utility 'gdxin' / '%input_dir%/ts_unit/' tSolve.tl:0 '.gdx';
44
45
        execute_load ts_unit_update=ts_unit;
        ts_unit(unit, *, f_solve(f), tt_forecast(t))
46
47
48
            ${  not mf_realization(mSolve, f) // Realization not updated
*                ts_unit_update(unit, *, f, t) // Update only existing values (zeroes need to be EPS)
                }
49
50
            = ts_unit_update(unit, *, f, t);
    ); // END if('ts_unit')
51

52
    // Update ts_effUnit
53
    if (mTimeseries_loop_read(mSolve, 'ts_effUnit'),
54
        put_utility 'gdxin' / '%input_dir%/ts_effUnit/' tSolve.tl:0 '.gdx';
55
        execute_load ts_effUnit_update=ts_effUnit;
56
57
58
59
60
        ts_effUnit(effGroupSelectorUnit(effSelector, unit, effSelector), *, f_solve(f), tt_forecast(t))
            ${  not mf_realization(mSolve, f) // Realization not updated
*                ts_effUnit_update(effSelector, unit, effSelector, *, f, t) // Update only existing values (zeroes need to be EPS)
                }
            = ts_effUnit_update(effSelector, unit, effSelector, *, f, t);
61
    ); // END if('ts_effUnit')
62
63

    // Update ts_effGroupUnit
64
    if (mTimeseries_loop_read(mSolve, 'ts_effGroupUnit'),
65
        put_utility 'gdxin' / '%input_dir%/ts_effGroupUnit/' tSolve.tl:0 '.gdx';
66
67
        execute_load ts_effGroupUnit_update=ts_effGroupUnit;
        ts_effGroupUnit(effSelector, unit, *, f_solve(f), tt_forecast(t))
68
69
70
            ${  not mf_realization(mSolve, f) // Realization not updated
*                ts_effGroupUnit_update(effSelector, unit, *, f, t) // Update only existing values (zeroes need to be EPS)
                }
71
72
73
            = ts_effGroupUnit_update(effSelector, unit, *, f, t);
    ); // END if('ts_effGroupUnit')
$offtext
74
75

    // Update ts_influx
76
    if (mTimeseries_loop_read(mSolve, 'ts_influx'),
77
        put_utility 'gdxin' / '%input_dir%/ts_influx/' tSolve.tl:0 '.gdx';
78
79
        execute_load ts_influx_update=ts_influx;
        ts_influx(gn(grid, node), f_solve(f), tt_forecast(t))
80
81
82
            ${  not mf_realization(mSolve, f) // Realization not updated
*                and ts_influx_update(grid, node, f, t) // Update only existing values (zeroes need to be EPS)
                }
83
84
            = ts_influx_update(grid, node, f, t);
    ); // END if('ts_influx')
85
86

    // Update ts_cf
87
    if (mTimeseries_loop_read(mSolve, 'ts_cf'),
88
        put_utility 'gdxin' / '%input_dir%/ts_cf/' tSolve.tl:0 '.gdx';
89
90
        execute_load ts_cf_update=ts_cf;
        ts_cf(flowNode(flow, node), f_solve(f), tt_forecast(t))
91
92
93
            ${  not mf_realization(mSolve, f) // Realization not updated
*                and ts_cf_update(flow, node, f, t) // Update only existing values (zeroes need to be EPS)
                }
94
95
            = ts_cf_update(flow, node, f, t);
    ); // END if('ts_cf')
96
97

    // Update ts_reserveDemand
98
    if (mTimeseries_loop_read(mSolve, 'ts_reserveDemand'),
99
        put_utility 'gdxin' / '%input_dir%/ts_reserveDemand/' tSolve.tl:0 '.gdx';
100
101
        execute_load ts_reserveDemand_update=ts_reserveDemand;
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f_solve(f), tt_forecast(t))
102
103
104
            ${  not mf_realization(mSolve, f) // Realization not updated
*                and ts_reserveDemand_update(restype, up_down, node, f, t) // Update only existing values (zeroes need to be EPS)
                }
105
106
            = ts_reserveDemand_update(restype, up_down, node, f, t);
    ); // END if('ts_reserveDemand')
107
108

    // Update ts_node
109
    if (mTimeseries_loop_read(mSolve, 'ts_node'),
110
        put_utility 'gdxin' / '%input_dir%/ts_node/' tSolve.tl:0 '.gdx';
111
112
        execute_load ts_node_update=ts_node;
        ts_node(gn(grid, node), param_gnBoundaryTypes, f_solve(f), tt_forecast(t))
113
114
115
            ${  not mf_realization(mSolve, f) // Realization not updated
*                and ts_node_update(grid, node, param_gnBoundaryTypes, f ,t) // Update only existing values (zeroes need to be EPS)
                }
116
117
118
119
120
            = ts_node_update(grid, node, param_gnBoundaryTypes, f, t);
    ); // END if('ts_node')

* --- NO FORECAST DIMENSION, SHOULD THESE BE HANDLED SEPARATELY? --------------
// Currently, only updated until the forecast horizon, but is this correct?
121
122

    // Update ts_fuelPriceChange
123
    if (mTimeseries_loop_read(mSolve, 'ts_fuelPriceChange'),
124
        put_utility 'gdxin' / '%input_dir%/ts_fuelPriceChange/' tSolve.tl:0 '.gdx';
125
126
        execute_load ts_fuelPriceChange_update=ts_fuelPriceChange;
        ts_fuelPriceChange(fuel, tt_forecast(t))
127
*            ${ ts_fuelPriceChange_update(fuel, t) } // Update only existing values (zeroes need to be EPS)
128
129
            = ts_fuelPriceChange_update(fuel, t);
    ); // END if('ts_fuelPriceChange')
130
131

    // Update ts_unavailability
132
    if (mTimeseries_loop_read(mSolve, 'ts_unavailability'),
133
        put_utility 'gdxin' / '%input_dir%/ts_unavailability/' tSolve.tl:0 '.gdx';
134
135
        execute_load ts_unavailability_update=ts_unavailability;
        ts_unavailability(unit, tt_forecast(t))
136
*            ${ ts_unavailability_update(unit, t) } // Update only existing values (zeroes need to be EPS)
137
138
            = ts_unavailability_update(unit, t);
    ); // END if('ts_unavailability')
139

140
    // Update the next forecast
141
    tForecastNext(mSolve)
142
        = tForecastNext(mSolve) + mSettings(mSolve, 't_forecastJump');
143
); // END if(tForecastNext)
144

145
* =============================================================================
146
* --- Optional forecast improvement -------------------------------------------
147
148
* =============================================================================

149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
if(mSettings(mSolve, 't_improveForecast'),
* Linear improvement of the central forecast towards the realized forecast,
* while preserving the difference between the central forecast and the
* remaining forecasts

    // Determine the set of improved time steps
    option clear = tt;
    tt(tt_forecast(t))
        ${ ord(t) <= tSolveFirst + mSettings(mSolve, 't_improveForecast') }
        = yes;

    // Temporary forecast displacement to reach the central forecast
    option clear = ddf;
    ddf(f_solve(f))
        ${ not mf_central(mSolve, f) }
        = sum(mf_central(mSolve, f_), ord(f_) - ord(f));

    // Temporary forecast displacement to reach the realized forecast
    option clear = ddf_;
    ddf_(f_solve(f))
        ${ not mf_realization(mSolve, f) }
        = sum(mf_realization(mSolve, f_), ord(f_) - ord(f));

* --- Calculate the other forecasts relative to the central one ---------------

    loop(f_solve(f)${ not mf_realization(mSolve, f) and not mf_central(mSolve, f) },
$ontext
* These don't work due to the wild cards in the parameter definitions!
        // ts_unit
        ts_unit(unit, *, f, tt(t))
            = ts_unit(unit, *, f, t) - ts_unit(unit, *, f+ddf(f), t);
        // ts_effUnit
        ts_effUnit(effGroupSelectorUnit(effSelector, unit, effSelector), *, f, tt(t))
            = ts_effUnit(effSelector, unit, effSelector, *, f, t) - ts_effUnit(effSelector, unit, effSelector, *, f+ddf(f), t);
        // ts_effGroupUnit
        ts_effGroupUnit(effSelector, unit, *, f, tt(t))
            = ts_effGroupUnit(effSelector, unit, *, f, t) - ts_effGroupUnit(effSelector, unit, *, f+ddf(f), t);
$offtext
        // ts_influx
        ts_influx(gn(grid, node), f, tt(t))
            = ts_influx(grid, node, f, t) - ts_influx(grid, node, f+ddf(f), t);
        // ts_cf
        ts_cf(flowNode(flow, node), f, tt(t))
            = ts_cf(flow, node, f, t) - ts_cf(flow, node, f+ddf(f), t);
        // ts_reserveDemand
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f, tt(t))
            = ts_reserveDemand(restype, up_down, node, f, t) - ts_reserveDemand(restype, up_down, node, f+ddf(f), t);
        // ts_node
        ts_node(gn(grid, node), param_gnBoundaryTypes, f, tt(t))
            = ts_node(grid, node, param_gnBoundaryTypes, f, t) - ts_node(grid, node, param_gnBoundaryTypes, f+ddf(f), t);
    ); // END loop(f_solve)

* --- Linear improvement of the central forecast ------------------------------

    loop(mf_central(mSolve, f),
$ontext
* These don't work due to the wild cards in the parameter definitions!
        // ts_unit
        ts_unit(unit, *, f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_unit(unit, *, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_unit(unit, *, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_effUnit
        ts_effUnit(effGroupSelectorUnit(effSelector, unit, effSelector), *, f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_effUnit(effSelector, unit, effSelector, *, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_effUnit(effSelector, unit, effSelector, *, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_effGroupUnit
        ts_effGroupUnit(effSelector, unit, *, f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_effGroupUnit(effSelector, unit, *, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_effGroupUnit(effSelector, unit, *, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
$offtext
        // ts_influx
        ts_influx(gn(grid, node), f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_influx(grid, node, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_influx(grid, node, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_cf
        ts_cf(flowNode(flow, node), f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_cf(flow, node, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_cf(flow, node, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_reserveDemand
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_reserveDemand(restype, up_down, node, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_reserveDemand(restype, up_down, node, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_node
        ts_node(gn(grid, node), param_gnBoundaryTypes, f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_node(grid, node, param_gnBoundaryTypes, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_node(grid, node, param_gnBoundaryTypes, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
    ); // END loop(mf_central)

* --- Recalculate the other forecasts based on the improved central forecast --

    loop(f_solve(f)${ not mf_realization(mSolve, f) and not mf_central(mSolve, f) },
$ontext
* These don't work due to the wild cards in the parameter definitions!
        // ts_unit
        ts_unit(unit, *, f, tt(t))
            = ts_unit(unit, *, f, t) + ts_unit(unit, *, f+ddf(f), t);
        // ts_effUnit
        ts_effUnit(effGroupSelectorUnit(effSelector, unit, effSelector), *, f, tt(t))
            = ts_effUnit(effSelector, unit, effSelector, *, f, t) + ts_effUnit(effSelector, unit, effSelector, *, f+ddf(f), t);
        // ts_effGroupUnit
        ts_effGroupUnit(effSelector, unit, *, f, tt(t))
            = ts_effGroupUnit(effSelector, unit, *, f, t) + ts_effGroupUnit(effSelector, unit, *, f+ddf(f), t);
$offtext
        // ts_influx
        ts_influx(gn(grid, node), f, tt(t))
            = ts_influx(grid, node, f, t) + ts_influx(grid, node, f+ddf(f), t);
        // ts_cf
        ts_cf(flowNode(flow, node), f, tt(t))
            = max(min(ts_cf(flow, node, f, t) + ts_cf(flow, node, f+ddf(f), t), 1), 0); // Ensure that capacity factor forecasts remain between 0-1
        // ts_reserveDemand
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f, tt(t))
            = max(ts_reserveDemand(restype, up_down, node, f, t) + ts_reserveDemand(restype, up_down, node, f+ddf(f), t), 0); // Ensure that reserve demand forecasts remains positive
        // ts_node
        ts_node(gn(grid, node), param_gnBoundaryTypes, f, tt(t))
            = ts_node(grid, node, param_gnBoundaryTypes, f, t) + ts_node(grid, node, param_gnBoundaryTypes, f+ddf(f), t);
    ); // END loop(f_solve)

); // END if(t_improveForecast)
288

289
290
291
292
293
294
295
296
297
298
299
300
301
302
* =============================================================================
* --- Aggregate time series data for the time intervals -----------------------
* =============================================================================

// Loop over the defined blocks of intervals
loop(cc(counter),

    // Retrieve interval block time steps
    option clear = tt_interval;
    tt_interval(t) = tt_block(counter, t);

    // If stepsPerInterval equals one, simply use all the steps within the block
    if(mInterval(mSolve, 'stepsPerInterval', counter) = 1,

303
304
* --- Select time series data matching the intervals, for stepsPerInterval = 1, this is trivial.

305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
        loop(ft(f_solve, tt_interval(t)),
            ts_cf_(flowNode(flow, node), f_solve, t, s)$msf(mSolve, s, f_solve)
                = ts_cf(flow, node, f_solve, t + (dt_sampleOffset(flow, node, 'ts_cf', s) + dt_circular(t)));
            ts_influx_(gn(grid, node), f_solve, t, s)$msf(mSolve, s, f_solve)
                = ts_influx(grid, node, f_solve, t + (dt_sampleOffset(grid, node, 'ts_influx', s) + dt_circular(t)));
            ts_unit_(unit, param_unit, f_solve, t)
              ${p_unit(unit, 'useTimeseries')} // Only include units that have timeseries attributed to them
                = ts_unit(unit, param_unit, f_solve, t+dt_circular(t));
            // Reserve demand relevant only up until reserve_length
            ts_reserveDemand_(restypeDirectionNode(restype, up_down, node), f_solve, t)
              ${ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')}
                = ts_reserveDemand(restype, up_down, node, f_solve, t+dt_circular(t));
            ts_node_(gn_state(grid, node), param_gnBoundaryTypes, f_solve, t, s)
              ${p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
                and msf(mSolve, s, f_solve)}
                = ts_node(grid, node, param_gnBoundaryTypes, f_solve, t + (dt_sampleOffset(grid, node, param_gnBoundaryTypes, s) + dt_circular(t)));
            // Fuel price time series
            ts_fuelPrice_(fuel, t)
                = ts_fuelPrice(fuel, t+dt_circular(t));
        ); // END loop(ft)

326
327
* --- If stepsPerInterval exceeds 1 (stepsPerInterval < 1 not defined) --------

328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
    elseif mInterval(mSolve, 'stepsPerInterval', counter) > 1,

        // Select and average time series data matching the intervals, for stepsPerInterval > 1
        // Loop over the t:s of the interval
        loop(ft(f_solve, tt_interval(t)),
            // Select t:s within the interval
            Option clear = tt;
            tt(t_)
                ${tt_interval(t_)
                  and ord(t_) >= ord(t)
                  and ord(t_) < ord(t) + mInterval(mSolve, 'stepsPerInterval', counter)
                 }
                = yes;
            ts_influx_(gn(grid, node), f_solve, t, s)$msf(mSolve, s, f_solve)
                = sum(tt(t_), ts_influx(grid, node, f_solve, t_ + (dt_sampleOffset(grid, node, 'ts_influx', s) + dt_circular(t_))))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ts_cf_(flowNode(flow, node), f_solve, t, s)$msf(mSolve, s, f_solve)
                = sum(tt(t_), ts_cf(flow, node, f_solve, t_ + (dt_sampleOffset(flow, node, 'ts_cf', s) + dt_circular(t_))))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ts_unit_(unit, param_unit, f_solve, t)
              ${ p_unit(unit, 'useTimeseries')} // Only include units with timeseries attributed to them
                = sum(tt(t_), ts_unit(unit, param_unit, f_solve, t_+dt_circular(t_)))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            // Reserves relevant only until reserve_length
            ts_reserveDemand_(restypeDirectionNode(restype, up_down, node), f_solve, t)
              ${ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')  }
                = sum(tt(t_), ts_reserveDemand(restype, up_down, node, f_solve, t_+dt_circular(t_)))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ts_node_(gn_state(grid, node), param_gnBoundaryTypes, f_solve, t, s)
              ${p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
                and msf(mSolve, s, f_solve)}
                   // Take average if not a limit type
                = (sum(tt(t_), ts_node(grid, node, param_gnBoundaryTypes, f_solve, t_ + (dt_sampleOffset(grid, node, param_gnBoundaryTypes, s) + dt_circular(t_))))
                    / mInterval(mSolve, 'stepsPerInterval', counter))$(not sameas(param_gnBoundaryTypes, 'upwardLimit') or sameas(param_gnBoundaryTypes, 'downwardLimit'))
                  // Maximum lower limit
                  + smax(tt(t_), ts_node(grid, node, param_gnBoundaryTypes, f_solve, t_ + (dt_sampleOffset(grid, node, param_gnBoundaryTypes, s) + dt_circular(t_))))
                      $sameas(param_gnBoundaryTypes, 'downwardLimit')
                  // Minimum upper limit
                  + smin(tt(t_), ts_node(grid, node, param_gnBoundaryTypes, f_solve, t_ + (dt_sampleOffset(grid, node, param_gnBoundaryTypes, s) + dt_circular(t_))))
                       $sameas(param_gnBoundaryTypes, 'upwardLimit');
            // Fuel price time series
            ts_fuelPrice_(fuel, t)
                = sum(tt(t_), ts_fuelPrice(fuel, t_+dt_circular(t_)))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ); // END loop(ft)

    ); // END if(stepsPerInterval)
); // END loop(counter)

Erkka Rinne's avatar
Erkka Rinne committed
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
* =============================================================================
* --- Input data processing ---------------------------------------------------
* =============================================================================

* --- Scenario reduction ------------------------------------------------------
s_active(s) = ms(mSolve, s);

if(active(mSolve, 'scenred'),
    $$include 'inc/scenred.gms'
);

* --- 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.

[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

* Influx
loop(gn(grid, node)$p_autocorrelation(grid, node, 'ts_influx'),
    ts_influx_mean(grid, node, ft(f, t))$mf_central(mSolve, f)
        = sum(s_parallel(s_active), ts_influx_(grid, node, f, t, s_active))
                / sum(s_parallel(s_active), 1);

    ts_influx_std(grid, node, ft(f, t))$mf_central(mSolve, f)
        = sqrt(sum(s_parallel(s_active), sqr(ts_influx_(grid, node, f, t, s_active)
                                         - ts_influx_mean(grid, node, f, t)))
                / sum(s_parallel(s_active), 1)
          );

    // Do smoothing
    loop(mst_end(ms_initial(mSolve, s_), t_),
        ts_influx_(grid, node, ft(f, t), s)$(ts_influx_std(grid, node, f, 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, f, t+dt_circular(t))
                    / ts_influx_std(grid, node, f, t_+dt_circular(t_)))
                * power(p_autocorrelation(grid, node, 'ts_influx'), abs(ord(t) - ord(t_)))
              ));
    );
);

* CF
loop(flowNode(flow, node)$p_autocorrelation(flow, node, 'ts_cf'),
    ts_cf_mean(flow, node, ft(f, t))$mf_central(mSolve, f)
        = sum(s_parallel(s_active), ts_cf_(flow, node, f, t, s_active))
                / sum(s_parallel(s_active), 1);

    ts_cf_std(flow, node, ft(f, t))$mf_central(mSolve, f)
        = sqrt(sum(s_parallel(s_active), sqr(ts_cf_(flow, node, f, t, s_active)
                                         - ts_cf_mean(flow, node, f, t)))
                / sum(s_parallel(s_active), 1)
          );

    // Do smoothing
    loop(mst_end(ms_initial(mSolve, s_), t_),
        ts_cf_(flow, node, ft(f, t), s)$(ts_cf_std(flow, node, f, 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, f, t+dt_circular(t))
                    / ts_cf_std(flow, node, f, t_+dt_circular(t_)))
                * power(p_autocorrelation(flow, node, 'ts_cf'), abs(ord(t) - ord(t_)))
              ));
    );
);