3c_inputsLoop.gms 27.3 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 (mSettings(mSolve, 'onlyExistingForecasts')
47
                   -> ts_unit_update(unit, param_unit, f, t)) // Update only existing values (zeroes need to be EPS)
48
                }
49
            = ts_unit_update(unit, param_unit, f, t);
50
    ); // END if('ts_unit')
51
52
53
54
$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
55
    // Update ts_effUnit
56
    if (mTimeseries_loop_read(mSolve, 'ts_effUnit'),
57
        put_utility 'gdxin' / '%input_dir%/ts_effUnit/' tSolve.tl:0 '.gdx';
58
        execute_load ts_effUnit_update=ts_effUnit;
59
        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
60
            ${  not mf_realization(mSolve, f) // Realization not updated
61
*               and ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t) // Update only existing values (zeroes need to be EPS)
62
                }
63
            = ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t);
64
    ); // END if('ts_effUnit')
65
66

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

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

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

    // Update ts_node
114
    if (mTimeseries_loop_read(mSolve, 'ts_node'),
115
        put_utility 'gdxin' / '%input_dir%/ts_node/' tSolve.tl:0 '.gdx';
116
117
        execute_load ts_node_update=ts_node;
        ts_node(gn(grid, node), param_gnBoundaryTypes, f_solve(f), tt_forecast(t))
118
            ${  not mf_realization(mSolve, f) // Realization not updated
119
                and (mSettings(mSolve, 'onlyExistingForecasts')
120
                     -> ts_node_update(grid, node, param_gnBoundaryTypes, f ,t)) // Update only existing values (zeroes need to be EPS)
121
                }
122
123
124
125
126
            = 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?
127
128

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

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

146
); // END if(tForecastNext)
147

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

292
293
294
295
296
297
298
299
300
301
* =============================================================================
* --- 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);
302
303
304
    // Make a temporary clone of tt_interval(t)
    option clear = tt_;
    tt_(tt_interval) = yes;
305
306
307
308

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

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

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

345
346
* --- If stepsPerInterval exceeds 1 (stepsPerInterval < 1 not defined) --------

347
348
349
350
    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
351
        loop(ft(f_solve(f), tt_interval(t)),
352
353
            // Select t:s within the interval
            Option clear = tt;
354
355
356
            tt(tt_(t_))
                ${  ord(t_) >= ord(t)
                    and ord(t_) < ord(t) + mInterval(mSolve, 'stepsPerInterval', counter)
357
358
                 }
                = yes;
359
360
361

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

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

423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
* --- 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
447
448
449
450
* =============================================================================
* --- Input data processing ---------------------------------------------------
* =============================================================================

451
$ifthen defined scenario
Erkka Rinne's avatar
Erkka Rinne committed
452
* --- Scenario reduction ------------------------------------------------------
453
if(active(mSolve, 'scenred') and mSettings('schedule', 'scenarios') > 1,
Erkka Rinne's avatar
Erkka Rinne committed
454
455
    $$include 'inc/scenred.gms'
);
456
$endif
Erkka Rinne's avatar
Erkka Rinne committed
457

458
459
460
461
462
463
464
465
* --- 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);


466
467
* --- Calculate sample displacements ------------------------------------------
Options clear = ds, clear = ds_state;
Erkka Rinne's avatar
Erkka Rinne committed
468
469
loop((mst_start(mSolve, s, t), ss(s, s_)),
    ds(s, t) = -(ord(s) - ord(s_));
470
471
472
473
474
475
    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
476
477
* --- Smooting of stochastic scenarios ----------------------------------------
$ontext
478
479
Smoothen the scenarios following the methodology presented in [1, p. 443].
This avoids a discontinuity `jump' after the initial sample.
Erkka Rinne's avatar
Erkka Rinne committed
480
481
482
483
484
485

[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

486
// Do smoothing
487
488
489
loop((ms_initial(mSolve, s_), t_, ft(f_, t__))
    $[ord(t_) = msEnd(mSolve, s_) and mst_end(mSolve, s_, t__)
      and (mf_realization(mSolve, f_) xor mf_central(mSolve, f_))],
490
491
    $$batinclude 'inc/smoothing.gms' ts_influx
    $$batinclude 'inc/smoothing.gms' ts_cf
Erkka Rinne's avatar
Erkka Rinne committed
492
);