3c_inputsLoop.gms 28.5 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
// Determine the necessary horizon for updating data
option clear = tmp;
28
29
tmp = max(  mSettings(mSolve, 't_forecastLengthUnchanging'),
            mSettings(mSolve, 't_forecastLengthDecreasesFrom')
30
31
32
33
34
35
36
            );

// 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) - mSettings(mSolve, 't_forecastJump'), // tForecastNext updated already in periodicLoop!
39

40
41
42
    // Update ts_unit
    if (mTimeseries_loop_read(mSolve, 'ts_unit'),
        put_utility 'gdxin' / '%input_dir%/ts_unit/' tSolve.tl:0 '.gdx';
43
        execute_load ts_unit_update=ts_unit;
44
        ts_unit(unit_timeseries(unit), param_unit, f_solve(f), tt_forecast(t)) // Only update if time series enabled for the unit
45
            ${  not mf_realization(mSolve, f) // Realization not updated
46
*               and ts_unit_update(unit, param_unit, f, t) // Update only existing values (zeroes need to be EPS)
47
                }
48
            = ts_unit_update(unit, param_unit, f, t);
49
    ); // END if('ts_unit')
50
51
52
53
$ontext
* !!! NOTE !!!
* These probably shouldn't be read at all, as p_effUnit and p_effGroupUnit are
* not input data, but calculated based on p_unit
54
    // Update ts_effUnit
55
    if (mTimeseries_loop_read(mSolve, 'ts_effUnit'),
56
        put_utility 'gdxin' / '%input_dir%/ts_effUnit/' tSolve.tl:0 '.gdx';
57
        execute_load ts_effUnit_update=ts_effUnit;
58
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f_solve(f), tt_forecast(t)) // Only update if time series enabled for the unit
59
            ${  not mf_realization(mSolve, f) // Realization not updated
60
*               and ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t) // Update only existing values (zeroes need to be EPS)
61
                }
62
            = ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t);
63
    ); // END if('ts_effUnit')
64
65

    // Update ts_effGroupUnit
66
    if (mTimeseries_loop_read(mSolve, 'ts_effGroupUnit'),
67
        put_utility 'gdxin' / '%input_dir%/ts_effGroupUnit/' tSolve.tl:0 '.gdx';
68
        execute_load ts_effGroupUnit_update=ts_effGroupUnit;
69
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f_solve(f), tt_forecast(t)) // Only update if time series enabled for the unit
70
            ${  not mf_realization(mSolve, f) // Realization not updated
71
*               and ts_effGroupUnit_update(effSelector, unit, param_eff, f, t) // Update only existing values (zeroes need to be EPS)
72
                }
73
            = ts_effGroupUnit_update(effSelector, unit, param_eff, f, t);
74
75
    ); // END if('ts_effGroupUnit')
$offtext
76
    // Update ts_influx
77
    if (mTimeseries_loop_read(mSolve, 'ts_influx'),
78
        put_utility 'gdxin' / '%input_dir%/ts_influx/' tSolve.tl:0 '.gdx';
79
80
        execute_load ts_influx_update=ts_influx;
        ts_influx(gn(grid, node), f_solve(f), tt_forecast(t))
81
82
83
            ${  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)
                }
84
85
            = ts_influx_update(grid, node, f, t);
    ); // END if('ts_influx')
86
87

    // Update ts_cf
88
    if (mTimeseries_loop_read(mSolve, 'ts_cf'),
89
        put_utility 'gdxin' / '%input_dir%/ts_cf/' tSolve.tl:0 '.gdx';
90
91
        execute_load ts_cf_update=ts_cf;
        ts_cf(flowNode(flow, node), f_solve(f), tt_forecast(t))
92
93
94
            ${  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)
                }
95
96
            = ts_cf_update(flow, node, f, t);
    ); // END if('ts_cf')
97
98

    // Update ts_reserveDemand
99
    if (mTimeseries_loop_read(mSolve, 'ts_reserveDemand'),
100
        put_utility 'gdxin' / '%input_dir%/ts_reserveDemand/' tSolve.tl:0 '.gdx';
101
102
        execute_load ts_reserveDemand_update=ts_reserveDemand;
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f_solve(f), tt_forecast(t))
103
104
105
            ${  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)
                }
106
107
            = ts_reserveDemand_update(restype, up_down, node, f, t);
    ); // END if('ts_reserveDemand')
108
109

    // Update ts_node
110
    if (mTimeseries_loop_read(mSolve, 'ts_node'),
111
        put_utility 'gdxin' / '%input_dir%/ts_node/' tSolve.tl:0 '.gdx';
112
113
        execute_load ts_node_update=ts_node;
        ts_node(gn(grid, node), param_gnBoundaryTypes, f_solve(f), tt_forecast(t))
114
115
116
            ${  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)
                }
117
118
119
120
121
            = 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?
122
123

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

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

141
); // END if(tForecastNext)
142

143
* =============================================================================
144
* --- Optional forecast improvement -------------------------------------------
145
146
* =============================================================================

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
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) },
        // ts_unit
174
175
176
177
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t))// Only update for units with time series enabled
            = ts_unit(unit, param_unit, f, t) - ts_unit(unit, param_unit, f+ddf(f), t);
$ontext
* Should these be handled here at all? See above note
178
        // ts_effUnit
179
180
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effUnit(effSelector, unit, effSelector, param_eff, f, t) - ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf(f), t);
181
        // ts_effGroupUnit
182
183
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effGroupUnit(effSelector, unit, param_eff, f, t) - ts_effGroupUnit(effSelector, unit, param_eff, f+ddf(f), t);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
$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),
        // ts_unit
203
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t)) // Only update for units with time series enabled
204
            = [ + (ord(t) - tSolveFirst)
205
                    * ts_unit(unit, param_unit, f, t)
206
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
207
                    * ts_unit(unit, param_unit, f+ddf_(f), t)
208
                ] / mSettings(mSolve, 't_improveForecast');
209
210
$ontext
* Should these be handled here at all? See above note
211
        // ts_effUnit
212
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
213
            = [ + (ord(t) - tSolveFirst)
214
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f, t)
215
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
216
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf_(f), t)
217
218
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_effGroupUnit
219
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
220
            = [ + (ord(t) - tSolveFirst)
221
                    * ts_effGroupUnit(effSelector, unit, param_eff, f, t)
222
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
223
                    * ts_effGroupUnit(effSelector, unit, param_eff, f+ddf_(f), t)
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
                ] / 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) },
        // ts_unit
260
261
262
263
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t)) // Only update for units with time series enabled
            = ts_unit(unit, param_unit, f, t) + ts_unit(unit, param_unit, f+ddf(f), t);
$ontext
* Should these be handled here at all? See above note
264
        // ts_effUnit
265
266
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effUnit(effSelector, unit, effSelector, param_eff, f, t) + ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf(f), t);
267
        // ts_effGroupUnit
268
269
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effGroupUnit(effSelector, unit, param_eff, f, t) + ts_effGroupUnit(effSelector, unit, param_eff, f+ddf(f), t);
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
$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)
286

287
288
289
290
291
292
293
294
295
296
* =============================================================================
* --- 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);
297
298
299
    // Make a temporary clone of tt_interval(t)
    option clear = tt_;
    tt_(tt_interval) = yes;
300
301
302
303

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

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

306
307
308
        loop(ft(f_solve(f), tt_interval(t)),
            ts_unit_(unit_timeseries(unit), param_unit, f, t) // Only if time series enabled for the unit
                = ts_unit(unit, param_unit, f, t+dt_circular(t));
309
310
311
312
313
314
315
$ontext
* Should these be handled here at all? See above comment
            ts_effUnit_(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f_solve, t)
                = ts_effUnit(effSelector, unit, effSelector, param_eff, f_solve, t+dt_circular(t));
            ts_effGroupUnit_(effSelector, unit_timeseries(unit), param_eff, f_solve, t)
                = ts_effGroupUnit(effSelector, unit, param_eff, f_solve, t+dt_circular(t));
$offtext
316
            ts_cf_(flowNode(flow, node), f, t, s)$sft(s, f, t)
317
318
                = ts_cf(flow, node, f + df_scenario(f, t), t + (dt_scenarioOffset(flow, node, 'ts_cf', s)
                                            + dt_circular(t)$(not gn_scenarios(flow, node, 'ts_cf'))));
319
            ts_influx_(gn(grid, node), f, t, s)$sft(s, f, t)
320
                = ts_influx(grid, node, f + df_scenario(f, t), t + (dt_scenarioOffset(grid, node, 'ts_influx', s)
321
                                                + dt_circular(t)$(not gn_scenarios(grid, node, 'ts_influx'))));
322
            // Reserve demand relevant only up until reserve_length
323
            ts_reserveDemand_(restypeDirectionNode(restype, up_down, node), f, t)
324
              ${ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')}
325
326
327
                = ts_reserveDemand(restype, up_down, node, f + df_scenario(f,t),
                                   t+dt_circular(t));
            ts_node_(gn_state(grid, node), param_gnBoundaryTypes, f, t, s)
328
              $p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
329
330
331
                = ts_node(grid, node, param_gnBoundaryTypes, f + df_scenario(f, t),
                          t + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                               + dt_circular(t)$(not gn_scenarios(grid, node, 'ts_node'))));
332
333
            // Fuel price time series
            ts_fuelPrice_(fuel, t)
334
              ${ p_fuelPrice(fuel, 'useTimeSeries') }
335
336
337
                = ts_fuelPrice(fuel, t+dt_circular(t));
        ); // END loop(ft)

338
339
* --- If stepsPerInterval exceeds 1 (stepsPerInterval < 1 not defined) --------

340
341
342
343
    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
344
        loop(ft(f_solve(f), tt_interval(t)),
345
346
            // Select t:s within the interval
            Option clear = tt;
347
348
349
            tt(tt_(t_))
                ${  ord(t_) >= ord(t)
                    and ord(t_) < ord(t) + mInterval(mSolve, 'stepsPerInterval', counter)
350
351
                 }
                = yes;
352
353
354

            ts_unit_(unit_timeseries(unit), param_unit, f, t)
                = sum(tt(t_), ts_unit(unit, param_unit, f, t_+dt_circular(t_)))
355
356
357
358
359
360
361
362
363
364
                    / 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, f_solve, t)
                = sum(tt(t_), ts_effUnit(effSelector, unit, effSelector, param_eff, f_solve, t_+dt_circular(t_))))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ts_effGroupUnit_(effSelector, unit_timeseries(unit), param_eff, f_solve, t)
                = sum(tt(t_), ts_effGroupUnit(effSelector, unit, param_eff, f_solve, t_+dt_circular(t_))))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
$offtext
365
366
367
368
            ts_influx_(gn(grid, node), f, t, s)$sft(s, f, t)
                = sum(tt(t_), ts_influx(grid, node, f + df_scenario(f, t),
                                        t_ + (dt_scenarioOffset(grid, node, 'ts_influx', s)
                                              + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_influx')))))
369
                    / mInterval(mSolve, 'stepsPerInterval', counter);
370
371
372
373
            ts_cf_(flowNode(flow, node), f, t, s)$sft(s, f, t)
                = sum(tt(t_), ts_cf(flow, node, f + df_scenario(f, t),
                                    t_ + (dt_scenarioOffset(flow, node, 'ts_cf', s)
                                          + dt_circular(t_)$(not gn_scenarios(flow, node, 'ts_cf')))))
374
375
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            // Reserves relevant only until reserve_length
376
            ts_reserveDemand_(restypeDirectionNode(restype, up_down, node), f, t)
377
              ${ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')  }
378
379
                = sum(tt(t_), ts_reserveDemand(restype, up_down, node,
                                               f + df_scenario(f, t), t_+dt_circular(t_)))
380
                    / mInterval(mSolve, 'stepsPerInterval', counter);
381
            ts_node_(gn_state(grid, node), param_gnBoundaryTypes, f, t, s)
382
              ${p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
383
                and sft(s, f, t)}
384
                   // Take average if not a limit type
385
386
387
                = (sum(tt(t_), ts_node(grid, node, param_gnBoundaryTypes,
                                       f + df_scenario(f, t), t_ + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                                                                    + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node')))))
388
                    / mInterval(mSolve, 'stepsPerInterval', counter))$(not (sameas(param_gnBoundaryTypes, 'upwardLimit') or sameas(param_gnBoundaryTypes, 'downwardLimit')))
389
                  // Maximum lower limit
390
391
392
                  + smax(tt(t_), ts_node(grid, node, param_gnBoundaryTypes,
                                         f + df_scenario(f, t), t_ + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                                                                      + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node')))))
393
394
                      $sameas(param_gnBoundaryTypes, 'downwardLimit')
                  // Minimum upper limit
395
396
397
398
                  + smin(tt(t_), ts_node(grid, node, param_gnBoundaryTypes,
                                         f + df_scenario(f, t),
                                         t_ + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                                               + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node')))))
399
400
401
                       $sameas(param_gnBoundaryTypes, 'upwardLimit');
            // Fuel price time series
            ts_fuelPrice_(fuel, t)
402
                ${ p_fuelPrice(fuel, 'useTimeSeries') }
403
404
405
406
407
408
409
                = sum(tt(t_), ts_fuelPrice(fuel, t_+dt_circular(t_)))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ); // END loop(ft)

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

410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
* --- Process unit time series data -------------------------------------------

// Calculate time series form parameters for units using direct input output conversion without online variable
// Always constant 'lb', 'rb', and 'section', so need only to define 'slope'.
loop(effGroupSelectorUnit(effDirectOff, unit_timeseries(unit), effDirectOff_),
    ts_effUnit(effDirectOff, unit, effDirectOff_, 'slope', ft(f, t))
        ${ sum(eff, ts_unit(unit, eff, f, t)) } // NOTE!!! Averages the slope over all available data.
        = sum(eff${ts_unit(unit, eff, f, t)}, 1 / ts_unit(unit, eff, f, t))
            / sum(eff${ts_unit(unit, eff, f, t)}, 1);
); // END loop(effGroupSelectorUnit)

// NOTE! Using the same methodology for the directOn and lambda approximations in time series form might require looping over ft(f,t) to find the min and max 'eff' and 'rb'
// Alternatively, one might require that the 'rb' is defined in a similar structure, so that the max 'rb' is located in the same index for all ft(f,t)

// Calculate unit wide parameters for each efficiency group
loop(effLevelGroupUnit(effLevel, effGroup, unit)${  mSettingsEff(mSolve, effLevel)
                                                    and p_unit(unit, 'useTimeseries')
                                                    },
    ts_effGroupUnit(effGroup, unit, 'lb', ft(f, t))${   sum(effSelector, ts_effUnit(effGroup, unit, effSelector, 'lb', f, t))}
        = smin(effSelector${effGroupSelectorUnit(effGroup, unit, effSelector)}, ts_effUnit(effGroup, unit, effSelector, 'lb', f, t));
    ts_effGroupUnit(effGroup, unit, 'slope', ft(f, t))${sum(effSelector, ts_effUnit(effGroup, unit, effSelector, 'slope', f, t))}
        = smin(effSelector$effGroupSelectorUnit(effGroup, unit, effSelector), ts_effUnit(effGroup, unit, effSelector, 'slope', f, t)); // Uses maximum efficiency for the group
); // END loop(effLevelGroupUnit)

Erkka Rinne's avatar
Erkka Rinne committed
434
435
436
437
438
439
440
441
442
* =============================================================================
* --- Input data processing ---------------------------------------------------
* =============================================================================

* --- Scenario reduction ------------------------------------------------------
if(active(mSolve, 'scenred'),
    $$include 'inc/scenred.gms'
);

443
444
445
446
447
448
449
450
* --- Update probabilities ----------------------------------------------------
Option clear = p_msft_probability;
p_msft_probability(msft(mSolve, s, f, t))
    = p_mfProbability(mSolve, f)
        / sum(f_${ft(f_, t)},
              p_mfProbability(mSolve, f_)) * p_msProbability(mSolve, s);


451
452
* --- Calculate sample displacements ------------------------------------------
Options clear = ds, clear = ds_state;
Erkka Rinne's avatar
Erkka Rinne committed
453
454
loop((mst_start(mSolve, s, t), ss(s, s_)),
    ds(s, t) = -(ord(s) - ord(s_));
455
456
457
458
459
460
    ds_state(gn_state(grid, node), s, t)
      ${not sum(s__, gnss_bound(grid, node, s__, s))
        and not sum(s__, gnss_bound(grid, node, s, s__))} = ds(s, t);
);


Erkka Rinne's avatar
Erkka Rinne committed
461
462
463
464
465
466
467
468
469
470
471
* --- 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

472
* Only do smoothing if scenarios are used
473
* TODO: Also enable when long-term forecast is used
474
if(mSettings(msolve, 'scenarios'),
475

Erkka Rinne's avatar
Erkka Rinne committed
476
477
* Influx
loop(gn(grid, node)$p_autocorrelation(grid, node, 'ts_influx'),
478
    ts_influx_mean(grid, node, t_active(t))$(sum(sft(s, f, t), 1))
479
480
        = sum(sft(s, f, t), ts_influx_(grid, node, f, t, s))
                / sum(sft(s, f, t), 1);
Erkka Rinne's avatar
Erkka Rinne committed
481

482
    ts_influx_std(grid, node, t_active(t))$(sum(sft(s, f, t), 1))
483
        = sqrt(sum(sft(s, f, t), sqr(ts_influx_(grid, node, f, t, s)
484
                                     - ts_influx_mean(grid, node, t)))
485
                / sum(sft(s, f, t), 1)
Erkka Rinne's avatar
Erkka Rinne committed
486
487
488
489
          );

    // Do smoothing
    loop(mst_end(ms_initial(mSolve, s_), t_),
490
        ts_influx_(grid, node, ft(f, t), s)$(ts_influx_std(grid, node, t_+dt_circular(t_))
Erkka Rinne's avatar
Erkka Rinne committed
491
492
493
494
495
496
                                             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))
497
498
                * (ts_influx_std(grid, node, t+dt_circular(t))
                    / ts_influx_std(grid, node, t_+dt_circular(t_)))
Erkka Rinne's avatar
Erkka Rinne committed
499
500
501
502
503
504
505
                * power(p_autocorrelation(grid, node, 'ts_influx'), abs(ord(t) - ord(t_)))
              ));
    );
);

* CF
loop(flowNode(flow, node)$p_autocorrelation(flow, node, 'ts_cf'),
506
    ts_cf_mean(flow, node, t_active(t))$(sum(sft(s, f, t), 1))
507
508
        = sum(sft(s, f, t), ts_cf_(flow, node, f, t, s))
                / sum(sft(s, f, t), 1);
Erkka Rinne's avatar
Erkka Rinne committed
509

510
    ts_cf_std(flow, node, t_active(t))$(sum(sft(s, f, t), 1))
511
        = sqrt(sum(sft(s, f, t), sqr(ts_cf_(flow, node, f, t, s)
512
                                     - ts_cf_mean(flow, node, t)))
513
                / sum(sft(s, f, t), 1)
Erkka Rinne's avatar
Erkka Rinne committed
514
515
516
          );

    // Do smoothing
517
518
    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_))
Erkka Rinne's avatar
Erkka Rinne committed
519
520
521
522
                                         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)
523
              + (ts_cf_(flow, node, f_, t_, s_)
Erkka Rinne's avatar
Erkka Rinne committed
524
                 - ts_cf_(flow, node, f, t_, s))
525
526
                * (ts_cf_std(flow, node, t+dt_circular(t))
                    / ts_cf_std(flow, node, t_+dt_circular(t_)))
Erkka Rinne's avatar
Erkka Rinne committed
527
528
529
530
                * power(p_autocorrelation(flow, node, 'ts_cf'), abs(ord(t) - ord(t_)))
              ));
    );
);
531

532
); // end if(mSettings(msolve, 'scenarios'),